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'An  analytical  and  experimental  study  has  been  carried  out  in  order  to  invest- 
igate the  behavior  of  a two-phase,  solid-liquid  system  ivs  which  the  solid  phase 
is  subjected  to  varying  heat  flux  and  temperature  at  its  boundaries. 

The  specific  geometry  chosen  for  the  study  was  that  of  a thin  solid  layer 
formation  on  a planar  surface  which  is  maintained  below  the  fusion  temperature 
of  the  liquid.  The  determination  of  the  phase  boundary  position  was  accomp- 
lished primarily  by  a variable  nodal  system  finite  difference  technique.  In 
addition,  an  integral  solution  similar  to  those  frequently  used  in  the  f' 
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20.  ABSTRACT  (Continued) 

^calculation  of  two-dimenaional  boundary  layer  flows  in  fluid  dynamics 
is  formulated  and  compared  with  the  results  of  the  finite  difference 
technique.  It  is  concluded  that  the  integral  solution  is  in  many 
instances  as  effective  a solution  technique  as  the  finite  difference 
method . 

An  experimental  study  was  carried  out  with  ice  and  water  in  turbulent 
channel  flow  in  order  to  test  the  validity  of  the  mathematical  model. 
The  results  of  the  experiment  indicate  that  the  model  is  capable  of 
predicting  the  position  of  the  phase  interface  to  within  limits  of 
accuracy  which  are  acceptable  for  most  engineering  applications.  It 
is  found  that  the  model's  inability  to  account  for  temperature 
dependent  solid  phase  properties  may  have  an  adverse  effect  on  its 
performance  for  substances  whose  transport  properties  are  highly 
dependent  upon  temperature,  ry 
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CHAPTER  I 


INTRODUCTION 

A . Discuss log  of  the  Problem  and  Mott vat  ion  f or  the  Work 

The  study  presented  here  was  carried  out  in  order  to  gain  a better 
understanding  of  the  phenomena  occurring  in  a two-phase,  solid-liquid 
system  subjected  to  varied  temperature  and  heat  flux  conditions  at  the 
system  boundaries.  In  this  case,  the  primary  objective  is  the  deter- 
mination of  the  phase  interface  position  as  a function  of  time  for  a 
given  set  of  boundary  and  initial  conditions. 

There  are  many  industrial  applications  in  which  some  knowledge 
of  the  behavior  of  such  a two-phase  system  may  be  a primary  factor  in 
the  design  or  operation  of  a given  apparatus.  Examples  of  such  applica- 
tions include  thermal  energy  storage,  food  processing,  metal  casting, 
water  desalinization,  and  the  design  and  operation  of  cryogenic  and 
liquid  metal  heat  exchangers.  The  solidification  problem  becomes 
especially  important  in  the  consideration  of  nuclear  reactor  safety 
where  molten  core  debris  may  solidify  and  block  coolant  flow  passages 
during  the  postaccidont  heat  removal  process  [see  Cheung  and  Baker  (1)]. 

There  are  several  factors  which  are  common  to  most  solidification 
and  molting  phenomena  of  practical  significance.  First,  there  is  a 
single,  continuous  phase  boundary  separating  the  solid  and  liquid  phases 
in  each  case.  The  shape  of  the  interface  (planar,  cylindrical,  etc.)  in 
general  depends  upon  the  specific  configuration  of  the  apparatus.  (There 
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are  some  situations  in  which  the  solid  phase  may  be  dispersed  throughout 
the  liquid  phase;  however,  this  situation  arises  only  when  the  liquid  is 
sufficiently  supercooled  and  will  not  be  considered  here.)  Also,  the 
boundary  conditions  on  the  solid  phase  are  such  that  there  is  one  fixed 
boundary  and  one  free  or  moving  boundary  where  the  moving  boundary  is  in 
direct  contact  with  the  liquid  phase  which  may  or  may  not  be  flowing  in 
forced  convection.  These  conditions  are  characteristic  of  the  class  of 
phase  change  phenomena  often  referred  to  as  "moving  boundary"  or 
"Stefan-type"  problems  and  describe  the  fundamental  physical  configuration 
to  be  studied  here. 


B.  Previous  Work  in  Phase  Change  Heat  Transfer 

The  basic  problem  of  heat  transfer  in  a substance  undergoing 
phase  change  has  avoided  general  solution  for  a long  time.  The 


analytical  difficulties  which  are  encountered  can  be  mainly  attributed 
to  the  nonlinear  character  of  the  interface  boundary  condition  associated 
with  the  governing  energy  equation.  Put  another  way,  the  spatial  domain 
on  which  the  solution  for  the  temperature  field  is  sought  changes 
continuously  with  time.  This  means  that  the  extent  of  each  phase  as 
well  as  the  temperature  distribution  in  the  phase  are  unknown  initially 
and  both  must  be  determined  simultaneously. 

The  techniques  of  solution  which  exist  at  this  time  are  basically 
of  three  general  types: 

Analytical  (Closed  Form)  Solutions  in  which  a particular  solution 
is  formulated  for  a set  of  simplified  boundary  conditions  (often 


rather  restrictive  and  unrealistic  from  a practical  point  of  view) 
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Approximate  Solutions  in  which  simplifying  assumptions  are  made, 
usually  in  order  to  reduce  the  number  of  independent  variables 
describing  the  system. 

Full  Finite  Difference  Techniques  in  which  the  energy  conservat ion 
equation  is  solved  in  its  most  general  form. 

A brief  description  of  each  method  is  given  below. 

1 . Analytical  Solutions 

Analytical  solutions  to  the  two-phase  conduction  problem  are 
available  for  only  a limited  number  of  cases  corresponding  to  very 
simple  boundary  conditions  and  arc  usually  of  limited  practical  interest. 
The  classical  example  of  a closed  form  solution  is  that  given  by  Stefan 
(2).  Stefan  considered  a semi-infinite  region  filled  with  the  isothermal 
(at  the  fusion  temperature)  liquid  phase,  bounded  by  a plane  wall  (fixed 
boundary).  At  zero  time,  the  wall  surface  temperature  is  stepped  down 
to  some  temperature  below  the  fusion  temperature  and  held  constant 
thereafter,  thus  initiating  solid  layer  formation  and  growth.  It  is 
assumed  that  heat  transfer  in  both  phases  is  via  conduction  only.  The 
solution  for  the  temperature  distribution  in  the  solid  is  obtained 
through  the  introduction  of  a similarity  variable  and  the  phase  boundary 
position  is  obtained  as  a simple,  algebraic  function  of  time. 

A somewhat  less  restrictive  result  was  obtained  by  Portnov  (3) 
for  the  instance  in  which  the  liquid  phase  is  not  initially  isothermal 
and  the  wall  temperature  is  a function  of  time.  The  solution  for  the 
semi-infinite  domain  is  obtained  through  Laplace  transform  methods  and 
ultimately  is  found  in  terms  of  an  infinite  series.  The  utility  of  this 


... 


solution,  however,  is  limited,  according  to  Stephan  (4),  due  to  the 
difficulties  encountered  in  the  evaluation  of  the  coefficients  in  the 
series. 
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It  is  apparent  from  just  these  two  examples  of  closed  form 
solutions  that  they  are  of  limited  applicability  and  that  other  methods 
must  be  developed  in  order  to  deal  with  problems  of  practical  importance. 

2 . Approximate  Techniques 

The  most  popular  approximate  techniques  for  solving  moving  boundary 
problems  are  the  integral  methods  in  which  the  temperature  profile  in 
each  phase  is  approximated  by  some  convenient  functional  relationship 
such  as  a polynomial  in  powers  of  the  spatial  coordinate.  The  energy 
equation  is  then  integrated  with  respect  to  the  space  variable  between 
the  moving  and  fixed  boundaries  using  the  approximate  profile.  This 
procedure  results  in  an  ordinary  differential  equation  (in  the  one- 
dimensional case)  for  the  phase  boundary  position  as  a function  of  time. 

The  procedure  described  here  parallels  the  integral  techniques  used  in 
the  analysis  of  two-dimensional  boundary  layer  flows  in  fluid  dynamics. 

Integral-type  analyses  have  been  applied  to  moving  boundary  systems 
by  Goodman  (5),  Libby  and  Chen  (6),  Savino  and  Siegel  (7),  and  Stephan  (4). 
The  rather  extensive  experimental  investigation  carried  out  by  Siegel 
and  Savino  yielded  results  which  were  in  good  agreement  with  their 
integral  formulation.  In  general,  the  integral  techniques  represent  a 
much  simplified  and  more  general  method  of  solution  in  comparison  to  the 
closed  form,  analytic  solutions  while  retaining  a high  degree  of  accuracy. 

Other  approximate  methods  which  have  been  employed  with  success 
in  phase  change  problems  are  the  variational  techniques  originally 
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developed  by  Biot  (8).  A variational  approach  was  applied  to  the  case 
of  constant  wall  temperature  solidification  by  Lapadula  and  Mueller  (9). 
This  work  resulted  in  solutions  for  the  phase  boundary  position  for  very 
small  and  very  large  time,  comparable  in  accuracy  to  the  integral  method 
used  by  Libby  and  Chen.  More  general  variational  techniques  (based  on 
the  Lagrangian  formalism)  were  employed  by  Bilenas  and  Jiji  (10)  in  a 
solution  of  the  free  boundary  solidification  problem  with  internal 
liquid  flow.  The  results  of  this  approximate  analysis,  however,  deviated 
as  much  as  +15%  from  the  supposedly  exact  full  finite  difference  solu- 
tions. 

On  the  whole,  variational  methods  are  not  as  popular  as  integral 
methods,  and  integral  methods  therefore  are  responsible  for  the  bulk  of 
the  approximate  solutions. 

3 . Full  Finite  Difference  Techniques 

The  third  major  category  of  solution  methods  are  the  numerical 
finite  difference  techniques.  The  principal  advantage  of  numerical 
solution,  of  course,  lies  in  the  capability  of  these  methods  to  handle 
irregular  geometries,  time  varying  boundary  conditions,  etc..  But  even 
finite  difference  methods  are  adversely  affected  by  the  moving  boundary 
characteristic  of  all  phase  change  phenomena.  In  fixed  nodal  arrange- 
ments, the  solid-liquid  interface  must  move  toward,  reach,  and  pass-by 
each  stationary  temperature  node  so  that  when  the  interface  is  located 
between  two  nodes,  its  exact  position  is  unknown  and  therefore  it  must 
be  located  approximately  by  interpolation.  A continuous  account  of 
fusion  front  travel  then  becomes  somewhat  of  a problem.  The  difference 
equations  and  interpolation  formulae  can  become  very  complex,  especially 


when  solution  in  two  space  variables  is  considered.  For  example,  see 
Springer  and  Olson  (11). 
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A unique  nodal  scheme  for  one-dimensional  problems  which  eliminates 
the  need  for  interpolation  and  yields  a continuous,  accurate  account  of 
the  fusion  front  motion  was  suggested  by  Murray  and  Landis  (12).  In 
this  nodal  arrangement,  the  nodal  mesh  changes  continuously  with  fusion 
front  travel  so  that  the  interface  position  is  always  coincident  with  a 
temperature  node;  thus,  there  is  no  need  for  interpolation  and  the 
difference  formulae  are  simplified  considerably. 

C.  The  Scope  and  Objectives  of  the  Present  Work 

The  purpose  of  the  present  study  is  twofold.  First,  it  is 
desired  to  formulate  a mathematical  model  which  may  be  used  to  predict 
the  phase  boundary  position  as  a function  of  time  in  a two-phase,  solid- 
liquid  system  similar  to  those  described  above.  The  model  accounts  for 
one-dimensional  heat  flow  in  the  solid  phase  and  completely  arbitrary  forms 
of  initial  and  boundary  conditions.  The  primary  method  of  analysis  is 
based  upon  a variable  nodal  system  finite  difference  method  in  the  solid 
phase  similar  to  that  used  by  Murray  and  Landis  (12).  In  the  formulation 
of  the  model,  it  is  assumed  that  the  mechanism  of  heat  transfer  to  the 
moving  phase  change  interface  is  by  convection  from  the  flowing  liquid 
because  this  condition  represents  a more  practically  important  case  than 
that  of  pure  conduction  heat  transfer  in  t ho  liquid.  In  addition,  an 
integral  solution  is  formulated  under  the  same  assumptions  as  those  made 
for  the  full  finite  difference  model.  The  two  solution  methods  are 
compared  both  qualitatively  and  quantitatively  in  order  to  determine  some 
of  the  advantages  and  disadvantages  of  integral  representation  as  opposed 


to  full  numerical  solution. 


Secondly,  an  experimental  investigat ion  is  carried  out  in  order 
to  test  the  analytical  procedures  discussed  above  for  boundary  and  initi- 
al conditions  representative  of  those  which  might  be  encountered  in 
actual  phase  change  phenomena.  The  experiment  incorporates  water  in 
turbulent  channel  flow  over  a planar  ice-water  interface.  The  transient 
growth  behavior  of  the  solid  phase  is  studied  and  compared  directly  with 
the  analytically  generated  solutions  for  the  phase  boundary  position. 

On  the  basis  of  this  comparison,  the  validity  and  applicability  of  the 
solution  techniques  is  evaluated. 


C 


CHAPTER  II 


ANALYSIS 


A.  Statement  of  the  Prob  lom 

Figure  1 shows  a two-phase  system  of  solid  and  liquid  in  which  the 
solid  layer  is  bounded  on  one  side  by  a planar,  rigid  wall  and  on  the 
other  side  by  the  flowing  liquid  phase.  Thu  wall  is  assumed  to  be 
either  a plane  wall  or  one  in  which  the  local  radius  of  curvature  is 
large  with  respect  to  the  local  solid  layer  thickness,  S.  At  the 
solid-liquid  interface,  y * S.  the  solid  phase  surface  is  subject  to  the 
convective  heat  flux  from  the  warmer  liquid  which  is  a function  of  the 
local  convective  conductance,  h,  and  the  freest  ream  temperature,  T . 

At  y ” 0,  the  layer  is  in  direct  contact  with  the  solid  wall  which  is 
assumed  always  to  be  at  some  temperature  lower  than  the  fusion  tempera- 
ture of  the  liquid;  hence,  some  solid  is  always  present.  The  assumpt ion 
of  direct  contact  between  the  solid  layer  and  wall  implies  the  absence 
of  any  contact  resistance  and,  therefore,  the  temperature  of  tho  solid 
at  y - 0 is  equal  to  the  wall  temperature.  It  is  important  to  note 
that,  in  this  analysis,  no  assumptions  will  be  made  regarding  the  form 
of  the  time  variations  of  wall  temperature  and  interior lal  heat  flux. 

On  the  other  hand,  both  T and  q”  must  be  known  as  t unctions  of  time 

w c 

before  the  problem  .solution  can  begin.  Specification  ot  the  heat  t lux 
at  the  interface  implies  that  the  fluid  dynamical  problem  tor  tho 
region  y N S is  Independent  of  the  solid  layer  thickness.  This  can 


q“ (HEAT  FLUX) 


(WALL) 


Figure  1.  Sketch  Illustrating  the  Basic  Arrangement 
of  Rigid  Wall,  Solid  Layer,  and  Liquid 
Freest  ream 
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usually  be  justified  as  a reasonable  assumption  in  external  flow  situa- 
tions when  S is  small,  but  for  internal  flows  when  the  layer  thickness 
becomes  appreciable  in  comparison  with  the  flow  radius,  this  assumption 
can  no  longer  be  made.  A more  complete  discussion  of  the  assumptions 
made  in  the  analysis  are  given  below. 

B.  Basic  Assumptions  in  Model  Formulation 

The  following  assumptions  are  made: 

1.  Conduction  heat  transfer  in  the  solid  layer  is  primarily  one- 
dimensional and  in  a direction  normal  to  the  solid  wall.  The  validity 
of  this  assumption  is  justified  for  thin  solid  layers  through  a compar- 
ison of  the  steady-state  heat  fluxes  in  the  y direction  and  in  the 
direction  parallel  to  the  solid  layer  surface  (near  the  wall).  This 
comparison  was  accomplished  experimentally  and  it  was  found  that  the 
flux  in  the  y direction  is  always  at  least  two  orders  of  magnitude 
greater  than  in  any  other  direction. 

2.  The  wall  temperature,  T , and  the  convective  heat  flux  at  the 

w 

interface,  q”,  are  not  dependent  upon  the  solid  layer  thickness. 

3.  The  material  under  consideration  is  a pure  substance.  The 
complexities  which  arise  due  to  phase  change  temperature  ranges  and 
solid  state  phase  transformations  in  multicomponent  systems  are  not 
considered  here.  With  respect  to  alloy  systems,  for  example,  Friedman 
(13)  has  approached  the  phase  change  problem  successfully  with  a finite 
element  technique;  however,  the  present  analysis  will  apply  only  to  pure 
substances,  homogeneous  in  the  solid  phase. 

A.  All  solid  properties  are  temperature  independent.  This 


assumption  should  not  severely  affect  the  accuracy  of  the  results  in 


problems  where  temperature  differences  are  not  extreme,  or  where 
properties  are  not  strongly  dependent  on  temperature. 
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5.  Thermodynamic  equilibrium  prevails  at  the  solid-liquid 


interface. 


C . Governing  Equations 


For  the  assumptions  stated  above,  conservation  of  energy  in  the 
solid  layer  and  at  the  phase  change  interface  yield  the  following  system: 


8 T _ 8T 
~ 2 ~ 8t  ’ 


0 < y < S 


T = To(y)  , 


S = S 

o 


t = 0 


T = T , 
w 


y = 0 


T = T , 
m 


y = s 


The  analytical  difficulties  which  arise  in  the  solution  of 
Equations  (1)  are  mainly  due  to  the  condition  (If)  at  the  phase  boundary. 
Consider  the  temperature  field  in  the  layer, 


T = T[y,t] 


si 
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At  the  interface i y - S, 


dT  - (~)dy  + (f~)dt  - 0 . 


Rearranging  and  substituting  into  condition  (If)  yields: 


. 3T  „ . .BT^T, 

k -r q “ -pit  l-r— /•»— ] 

Oy  nc  l3t  3y 


The  interface  boundary  condition  is  therefore  seen  to  be  nonlinear. 

The  system  (1)  is  nondimensional ised  by  defining  the  following 
dimens  ion  less  groups : 


T-T 

. y .tot  - w ...  s 

y*  s • 1 72  ■ 0 " rr  • s ' r • 

o S m w o 

o 


where  S >0  denotes  the  solid  layer  thickness  at  zero  time, 
o 

t - t*  - 0 . 

Substituting  into  Equations  (1)  results  in, 


_a_L  „ _L  dlt.  ,0_n  + jo. 

3y*2  st  dt*  1 1 at*  ’ 


0 < y*  < S* 


0 - 0o(y*)  , 


S*  - 1 , 


t*  - 0 


0-0  , 


y*  - 0 


— 


I 1 


and 


0-1 


30  , 

St  9y*  - 


dS* 

dt'*  ’ 


y*  - S* 


(V) 


(■>1) 


where  the  parameter,  St  (Stefan  number),  is  defined  by 


St 


e(T  -T  ) 
m w 


(6) 


and  tin'  dimensionless  heat  flux  parameter  by: 


<i* 


S 

o 

P C» 


(7) 


Equations  (5)  completely  define  the  problem  studied  here.  It  has 
been  shown  that  the  condition  (5f)  renders  the  system  nonlinear  and 
thereby  makes  closed  form  analytical  solutions  very  difficult  to  obtain 
for  all  but  the  most  elementary  forms  of  boundary  and  Initial  conditions. 

The  following  sections  briefly  describe  three  methods  of  solution,  of 
varying  degrees  of  generality,  which  under  the  appropriate  circumstaj\ces  # . 

are  found  to  yield  accurate  results  for  the  phase  change  interface 
position,  S*(t*)  . 

The  first  solution  (asymptotic)  is  a result  of  neglecting  the 
specific  heat  capacity  of  the  solid  phase  in  comparison  to  the  latent 
heat  of  fusion.  This  solution,  valid  only  for  small  values  of  the 
parameter,  St,  is  available  in  closed  form.  The  second  method  of 
solution  (Integral)  is  approximate,  but  more  general  than  the  asymptotic 
solution.  In  the  integral  technique,  0 and  y*  are  essentially 
eliminated  from  Equation  (5a) . and  the  resulting  ordinary  differential 


m llnlll 
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equation  for  S*  as  a function  of  time  is  readily  solved.  The  third 
solution  method  involves  the  finite  difference  representation  of 
Equations  (5)  and  is  completely  general  with  respect  to  the  form  of  the 
boundary  conditions  at  y*  = 0 and  y*  = S*  . 

D.  Asymptotic  Solution  for  Small  Stefan  Number 

The  parameter,  St  (Stefan  number)  which  arises  from  the  non- 
dimensional  ization  of  Equation  (5a)  is  an  important  parameter  in 
determining  the  way  heat  is  transferred  within  the  solid  layer  during 
layer  growth.  Qualitatively,  the  Stefan  number  is  often  interpreted  as 
a ratio  of  sensible  heat  capacity  effects  to  latent  heat  effects  in  the 
layer.  If  the  Stefan  number  is  small,  the  rate  of  the  solidification  or 
melting  process  is  essentially  controlled  by  the  latent  heat  of  fusion, 
i . In  this  instance,  the  temperature  profile  in  the  layer  will  respond 
more  quickly  to  a change  in  wall  temperature  than  will  the  layer  thickness 
and  the  profile  will  assume  a nearly  linear  shape  throughout  the  transient. 
If  the  assumption  of  a quasisteady  or  linear  temperature  profile  in  the 
layer  is  made,  a closed  form,  analytic  solution  for  the  phase  boundary 
position  becomes  possible. 

With  the  linear  assumption, 


0 


(8) 


Using  this  result  in  the  interface  heat  balance  expression,  (5f),  results 
in 


- q*  . 


(9) 


dS*  _ St 
dt*  S* 


If  it  is  assumed  further  that  wall  temperature  and  interface  heat  flux 
renvi in  constant,  then  rearrangement  und  integration  of  Equation  (9) 
using  condition  (5c)  yields 


t* 


St-q*t) 


(10) 


Equation  (10)  is  an  expression  for  the  phase  boundary  position  valid  for 
the  case  of  small  Stefan  number  (small  heat  capacity).  The  extent  to 
which  Equation  (10)  yields  accurate  results  may  only  be  established  by 
direct  comparison  with  an  exact  solution  for  S*(t*)  . 


E . Approximate  Solution  by  an  Integral  Technique 

The  system,  (5),  may  be  solved  approximately  by  estimating  the 
temperature  profile  in  the  solid  layer  with  some  convenient  algebraic 
expression  which  satisfies  the  boundary  conditions  exactly  at  y “ 0 
and  y a S and  the  differential  equation,  (5a),  only  in  the  mean  [see 
(14)1  f°r  the  region  0 < y < S . Equation  (5a)  is  thereby  reduced  to 
an  ordinary  differential  equation  for  S,  the  solid  layer  thickness,  as 
a function  of  time  which  may  be  easily  solved  by  standard  techniques, 
integrating  Equation  (5a)  over  the  solid  layer  yields: 


4 
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results  in: 


dS*  . * 30 

+ o*  - 


dt* 


dt* 


y*=0 


[' 


S* 

St ]|  0dy*-S* 


(13) 


which  is  the  basic  form  of  the  integral  equation  for  the  solid  layer. 

The  form  of  the  algebraic  expression  chosen  to  represent  the  solid 
layer  temperature  profile  is  a second  degree  polynomial  in  the  dimension- 
less position  coordinate  y*/S*  : 


6 ‘ Vsl>2  + Vs?  * D3  • 

The  constants  D^,  D2>  and  are  chosen  so  that  the  temperature 

profile  satisfies  the  boundary  conditions: 


(14) 


0 = 0 


y * = 0 


(15) 


0 = 1 


_ 30  . dS* 

St  “ q dT* 


Then  Equation  (14)  becomes: 


y*  = S* 


[*  - £ m ♦ <*>]  <s> * [f  c + "*> 


(16) 


(17) 


-i]  (r«)2  •(18) 


Substituting  Equation  (18)  into  Equation  (13)  results  in: 


6(sl  _ _ q*>  = -A-  sts*  + (~J  + q*) 

VS*  dt*  H ' dt*  2 dt* 


dt* 


[“ 


•] 


(19a) 
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subject  to  the  Initial  conditions 


and 


S* 


dS* 

dt* 


St 


90 

9y* 


y*  = S* 


t*=0 


(19b) 


(19c) 


The  system,  (19),  was  solved  numerically  via  a fifth-order  Runge- 
Kutta  numerical  integration  formula  utilizing  Sarafyan  embedding  for  step 
size  control  (15). 

F.  Finite  Difference  Solution 

The  finite  difference  solution  of  Equations  (5)  is  achieved  by 
employing  variable  nodal  spacing  in  the  solid  layer  (Figure  2)  similar 
to  that  used  by  Murray  and  Landis  (12).  This  procedure  greatly  simplifies 
the  solution  in  the  vicinity  of  the  moving  phase  boundary  because  the 
position  of  the  interface  always  coincides  with  one  of  the  temperature 
nodes  and  interpolation  for  the  phase  boundary  position  becomes 
unnecessary.  In  order  to  develop  an  equation  applicable  to  the  moving 
nodal  point  system,  the  time  derivative  in  (5a),  90/9t*,  is  expressed 

in  terms  of  the  total  time  derivative,  D0/Dt*,  and  the  velocity  of  the 
moving  nodes,  ^ — - . Thus, 


90  = _D0_  2*  d§*  90 

9t*  - Dt*  - S*  dt*  9y* 


(20) 


Substituting  Equation  (20)  in  Equation  (5a)  yields. 


920 

9y*2 


_ rn.n  + HL  d£i  _90_ 

St  dt*  1 J S*  dt*  9y* 


DO 

Dt* 


0 < y*  < S*  (21a) 
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with  the  associated  conditions. 


0-0  (y *) 

o 


t*  - 0 


(21b) 


S*  - 1 


(21c) 


0-0 


y*  - 0 


( 2 1 d ) 


0 = 1 


(2  le) 


...  90  . dS* 

3y*  ^ dt* 


( 2 1 f ) 


Equation  (21a)  is  solved  for  the  temperature  distribution  in  the 
solid  layer  by  an  explicit  forward  difference  technique.  The  resulting 
interfacial  temperature  gradient  30/3y*  j is  approximated  via  a 

three-point  backward  difference  formula  and  used  in  condition  (21f)  to 
yield  the  phase  boundary  velocity  and,  ultimately,  the  interface 
position.  The  detailed  difference  equations  can  be  found  in  Appendix  A. 

The  explicit  form  of  the  difference  equations  simplifies  the 
numerical  solution  of  Equation  (21a)  considerably;  however,  with  an 
explicit  solution  scheme,  the  question  of  stability  must  be  considered. 

II 

According  to  Ozlsik  (16),  there  is  no  general  stability  criterion  for 
determining  the  stability  limits  in  nonlinear  problems.  In  most  cases. 


the  stability  bounds  are  determined  by  numerical  experiment.  In  the 


present  study,  both  the  accuracy  and  stability  of  the  solution  are 
monitored  bv  continuously  correcting  the  temperature  distribution  and 
the  solid  layer  thickness  at  each  time  step.  The  corrections  are 
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applied  until  successive  corrections  agree  to  within  some  specified 
tolerance.  If  this  criterion  is  not  met  after  four  corrections,  the 
time  step  is  halved  and  that  time  stop  repeated.  In  all  the  runs  made 
(over  one  hundred) , there  were  no  instances  of  instability. 

The  results  of  the  asymptotic  solution  for  small  St,  the  integral 
technique,  and  the  finite  difference  model  are  illustrated  and  discussed 
fully  in  Chapter  IV. 

[ 


CHAPTER  III 


EXPERIMENTAL  INVESTIGATION 

A.  Introduction  and  Purpose  of  the  Experiment 

The  major  objective  in  carrying  out  this  experimental  investigation 
was  to  demonstrate  the  validity  and  the  range  of  applicability  of  the 
general  analytical  model  and  numerical  solution  procedure  discussed  in 
Chapter  II.  The  experiment  should  provide,  above  all,  a means  of 
evaluating  the  generality  of  the  theoretical  model  with  respect  to 
arbitrary  time  varying  conditions  at  the  solid  layer  boundaries.  On  the 
other  hand,  the  assumption  of  uni-dimensional  heat  flow  in  the  solid- 
layer,  and  the  assumptions  that  the  material  in  question  is  a pure 
substance  and  that  thermodynamic  equilibrium  exists  at  the  phase  change 
interface  must  be  preserved  in  the  experiment  in  order  to  honor  the 
assumptions  under  which  the  model  was  formulated.  In  particular,  the 
assumption  of  one-dimensionality  is  essential  to  the  numerical  procedure 
used  in  Chapter  II  to  solve  the  governing  energy  equation.  The  variable 
space  node  method  adopted  for  this  purpose  becomes  extremely  complex  if 
multidimensional  heat  flow  is  considered.  To  the  author's  knowledge,  the 
variable  nodal  system  has  not  been  implemented  in  multidimensional  phase 
change  problems  to  date.  Thus,  the  experimental  apparatus  used  in  this 
investigation  was  designed  in  such  a way  as  to  provide  a means  of 
evaluating  the  theoretical  model  without  violating  the  basic  assumptions 
under  which  the  problem  was  formulated. 
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In  order  to  effect  the  comparison  between  experiment  and  theory, 
several  basic  physical  measurements  must  be  made  with  a reasonable  degree 
of  accuracy  in  the  experiment.  The  most  important  measurements  to  be 
made  are  suggested  by  Figure  1 and  include  the  local  convective  heat 
flux  at  the  solid-liquid  interface  (y  = S) , the  local  temperature  at  the 
interface  between  cooling  surface  and  solid  layer  (y  ■ 0),  and  the  local 
solid  layer  thickness,  S . The  first  two  measurements  are  the  independ- 
ent or  driving  variables  of  the  solid  layer  growth  process,  while  the 
layer  thickness  is  the  resultant  of  these.  Therefore,  no  effort  was 
expended  in  the  experimental  design  to  provide  precise  control  over  the 

wall  temperature,  T , and  heat  flux  at  the  interlace,  q" , as  long  as 

W c 

the  variation  of  each  with  time  could  be  accurately  monitored  and 
recorded.  It  then  becomes  a simple  nutter  to  use  these  recorded 
variations  in  the  analytical  model  and  thus  generate  theoretical  solid 
thickness  variations  which  may  be  compared  with  the  experimental  results. 
On  the  basis  of  this  comparison,  the  solution  procedures  set  forth  in 
Chapter  II  can  be  evaluated. 

B . Experimental  Apparatus 

1.  Requirement  s 

The  experimental  apparatus  must,  to  a reasonable  degree,  duplicate 
the  conditions  shown  in  Figure  1.  Therefore,  the  geometry  chosen  for  the 
experiment  was  that  of  a plane  surface  over  which  the  liquid  phase  could 
flow  in  forced  convection.  The  flow  regime  of  the  liquid  flow  past  the 
solid  layer  was  chosen  to  be  one  of  parallel,  turbulent  flow.  The 
distinction  made  here  is  important  because,  to  the  author's  knowledge. 


there  has  been  no  other  experimental  study  nude  of  a moving  boundary 


23 


system  which  employed  turbulent  flow  .it  the  phase  change  interface. 

Also,  precautions  were  taken  which  would  cause  the  ice  layer  to  remain 
relatively  thin  so  that  heat  flow  within  it  would  remain  essentially 
one-dimensional.  The  presumption  that  the  heat  conduction  in  the  layer 
is  or.e-dimensional  would  be  verified  experimentally.  Provisions  also 
were  made  to  provide  a range  of  heat  flux  conditions  at  the  interface. 

This  was  accomplished  by  varied  liquid  temperature  and/or  liquid  flow 
velocity.  The  determination  of  the  interface  heat  flux  can  be  approached 
from  two  basic  directions.  Theories  of  convection  heat  transfer  may  be 
applied  to  predict  the  local  film  coefficient  at  a particular  point  on 
the  solid  layer  surface,  or  the  heat  flux  rany  be  determined  experimentally 
using  the  steady-state  solid  layer  thickness  as  a form  of  "heat  flux 
meter."  The  latter  approach  is  the  one  adopted  here,  while  convection 
theory  is  used  to  more  or  less  qualitatively  verify  the  results.  This 
approach  was  taken  because  the  convection  heat  transfer  situation  in  the 
experimental  apparatus  does  not  correspond  very  closely  to  any  well-tested 
heat  transfer  solution  in  the  literature.  Therefore,  the  experimental 
approach  should  prove  to  yield  much  more  reliable  estimates  of  the  inter- 
face heat  flux  and  thus  provide  better  information  for  the  analytical 
model.  In  addition  to  a variable  interface  heat  flux,  the  experimental 
apparatus  must  provide  a means  of  varying  the  solid  wall  temperature 
(but  not  necessarily  controlling  it).  All  the  primary  measurements 
mentioned  in  this  section  must  be  made  dynamically  if  solid  layer  growth 


transients  are  to  bo  satisfactorily  monitored. 
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2.  Basic  Configuration 

In  order  to  effect  the  above  requirements,  a closed  flow  loop  was 
constructed.  The  working  fluid  chosen  to  be  used  in  the  flow  loop  was 
ordinary  tap  water  because  of  its  high  availability  and  convenient 
physical  properties  (transparency,  moderate  fusion  temperature,  etc.). 

Other  fluids  were  considered,  including  low  melting  point  eutectic  alloys 
composed  of  bismuth,  lead,  tin,  and  cadmium.  However,  these  materials 
presented  rather  substantial  difficulties  in  terras  of  experimental 
measurements.  The  greatest  problem  being  their  opaqueness  which  would 
nuke  solid  layer  thickness  measurements  very  difficult.  Other  disadvantages 
of  the  low  melting  point  alloys  are  their  high  density , above  ambient 
fusion  temperatures,  and  the  extremely  high  fluidity  of  these  substances 
in  the  liquid  phase. 

The  basic  components  of  the  flow  loop  are  shown  in  Figure  3. 

These  consist  of  various  tanks  for  fluid  temperature  control,  pumps  for 
forc.d  circulation  and  cooling,  a rectangular  plexiglas  flow  channel 
incorporating  the  test  section  and  cooling  plate,  a flow  metering  device, 
and  the  interconnecting  tubing.  All  tubing  was  made  from  2.54  cm  (1  inch) 
inside  diameter  water  hose.  The  relatively  large  diameter  was  chosen  in 
order  to  minimize  pressure  loss  through  the  system  resulting  from  the 
high  velocities  and  volumetric  flow  rates  necessary  to  insure  turbulent 
flow  in  the  test  section.  Flow  Reynolds  numbers  used  in  the  experiment 
ranged  as  high  as  24,000  based  upon  the  flow  channel  hydraulic  diameter. 

The  cooling  plate,  the  leading  edge  of  which  was  a distance  of  64  cm 
(25.3  inches)  downstream  of  the  flow  channel  entrance,  was  oriented 
horizontally  and  formed  the  upper  wall  of  the  test  section.  Cooling 
requirements  for  the  plate  were  found  to  be  substantial  in  terms  of 
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temperature  and  total  heat  load  due  to  the  very  efficient  energy  trans- 
port characteristics  of  the  turbulent  channel  flow.  Cooling  of  the 
plate  was  accomplished  totally  by  means  of  a dry  ice-methyl  alcohol  bath 
located  immediately  above  the  cooling  plate.  A detailed  description  of 
each  component  in  the  loop  and  its  particular  function  is  given  in  the 
following  sections. 

3 . Flow  Channe l an d Cool ing  P late 

The  primary  piece  of  apparatus  in  the  system  was  the  clear  plexi- 

glas  flow  channel  housing  the  test  section  and  cooling  plate,  see  Figure 

4.  The  flow  channel  was  fabricated  from  a U-shaped  block  of  clear 

plexiglas  92.7  cm  (37.5  inches)  in  length  with  an  outside  width  of  7.62  cm 

(3  inches).  The  channel  was  sealed  with  a clear  plexiglas  cover  which 

was  through-bolted  to  the  channel  with  18  bolts,  each  12.7  cm  (5  inches) 

in  length.  The  interface  between  the  U-channel  and  the  cover  was  sealed 

with  a multipiece  cork  gasket  of  0.32  cm  (1/8  inch)  thickness.  The 

channel  Inside  dimensions  made  it  nearly  square  (3.75  cm  high  by  3. SI  cm 

2 

wide)  with  a total  flow  area  of  14.29  cm  . Tills  sizeable  flow  area 
coupled  with  the  requirement  of  turbulent  flow  past  the  cooling  plate 
resulted  in  substantial  mass  flow  rates  through  the  system.  The  channel 
inlet  had  a geometry  similar  to  an  abrupt  expansion  which  assured  the 
existence  of  a turbulent  boundary  layer  at  the  plate  leading  edge  even  in 
the  low  Reynolds  number  flow  ranges.  The  cooling  plate  was  fabricated 
from  pure  copper  and  was  a uniform  0.508  cm  (0.20  inches)  in  thickness. 

The  pi. iti'  was  18.8  cm  (7.4  inches)  long  in  the  direction  of  flow  and 
spanned  the  entire  channel  wfdili  forming  the  upper  wall  of  the  flow 
channel.  Copper  was  chosen  as  the  plate  material  for  several  reasons. 


First,  the  high  thermal  conductivity  of  copper  would  serve  to  minimize 
temperature  drop  through  the  plate  for  a given  convective  heat  flux  at 
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the  solid-liquid  interface.  Small  temperature  gradients  through  the 
plate  were  desirable  in  order  to  minimize  errors  in  plate  surface 
temperature  measurement,  this  point  will  be  discussed  in  greater  detail 
in  a later  section.  Also,  the  high  conductivity  and  thermal  diffusivity 
of  copper  served  to  minimize  the  temperature  gradients  in  directions 
transverse  to  that  of  ice  layer  growth.  Thus,  the  desired  one-dimension- 
ality of  the  heat  conduction  in  the  ice  layer  would  be  very  nearly 
achieved.  A further  check  on  the  one-dimensional  assumption  was  made 
experimentally  (see  "Steady-State  Heat  Flux  Measurements").  The  position 
of  the  cooling  plate  leading  edge  with  respect  to  the  flow  channel  inlet 
is  such  that  t he  turbulent  flow  is  for  all  practical  purposes  hydro- 
dynamieally  fully  developed  where  cooling  begins;  however,  at  the 
position  on  the  plate  at  which  ice  thickness  measurements  were  made,  the 
boundary  layer  was  somewhat  underdeveloped  thermally  and  this  fact  is 
confirmed  by  the  steady-state  heat  flux  measurements.  As  was  mentioned 
previously,  plate  cooling  was  achieved  by  means  of  a dry  ice-methyl 
alcohol  bath  which  filled  a bin  above  the  cooling  plate.  The  bin  was 
fabricated  from  brass  plate  which  was  then  silver  soldered  directly  to 
the  upper  surface  of  the  cooling  plate.  The  hath  had  an  approximate 
capacity  of  2.5  liters  which  was  found  to  he  sufficient  to  meet  the 
cooling  requirements  of  the  average  test  run.  The  bath  required  no 
external  agitation  as  sufficient  mixing  was  achieved  through  the  sub- 
limation process  of  the  dry  ice  itself.  On  the  other  hand,  it  was 
found  necessary  to  periodically  change  the  alcohol  in  the  bin  because 
it  would  absorb  moisture  from  the  surroundings  and,  over  a period  of 
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time,  the  water  dissolved  in  the  alcohol  would  affect  the  ability  of  the 
bath  mixture  to  reach  sufficiently  low  temperatures.  The  plate  surface 
temperature  (ice  layer  side)  produced  by  the  bath  was  in  general  a 
function  of  the  rate  of  heat  flow  through  the  cooling  plate,  which  in 
turn  was  a function  of  the  ice-water  interface  heat  flux  under  steady- 
state  ice  layer  conditions.  Visual  observation  of  the  ice  layer  surface 
and  the  layer  growth  process  was  enhanced  by  fine  polishing  the  vertical 
walls  of  the  flow  channel  in  the  test  section  area.  Ice  layer  observation 
was  also  improved  by  back-lighting  the  test  section  at  an  angle  which 
caused  light  to  reflect  off  the  ice  layer  surface.  The  entire  flow 
channel,  excluding  the  observation  windows  in  the  vicinity  of  the  test 
section,  was  thermally  insulated  with  3. SI  cm  (1.5  inch)  thick  poly- 
urethane foam  insulation.  This  heavy  insulation  assured  that  the  flow 
upstream  of  the  cooling  plate  would  be  uniform  in  temperature.  Thus, 
channel  centerline  temperatures  measured  upstream  of  the  test  section 
would  be  equivalent  to  mixed  mean  fluid  temperature.  As  the  water  passed 
through  the  test  section  and  passed  the  ice  layer,  the  fluid  bulk  or 
mixed  mean  temperature  was  considered  to  remain  essentially  constant  and 
equal  to  the  inlet  temperature  due  to  the  relatively  small  amount  of 
heat  transferred  to  the  water  per  unit  mass  flowing.  This  assumption  is 
verified  experimentally  by  measuring  the  bulk  water  temperature  down- 
stream of  the  test  section.  For  this  purpose,  the  flow  channel  and 
connecting  tubing  downstream  of  the  test  section  was  also  heavily 
insulated.  The  entire  flow  channel  assembly  was  rigidly  bolted  to  a 
steel  frame  so  that  the  test  section  observation  windows  were  supported 
at  an  elevation  which  provided  for  easy  visual  inspection  of  the  ice-water 
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interface.  Figure  5 is  a photograph  of  the  entire  flow  channel  assembly 
with  the  insulation  partially  removed. 

4.  Mixing  Tank 

A cylindrical  stainless  steel  mixing  tank  was  located  just  down- 
stream of  the  flow  channel  exit  and  test  section.  The  tank,  measuring 
17.8  cm  (7  inches)  in  diameter  and  25.4  cm  (10  inches)  in  length,  was 
baffled  so  that  the  water  flowing  through  it  would  be  completely  mixed 
when  it  reached  the  tank  exit  where  the  water  temperature  was  measured 
by  a thermocouple  probe.  In  this  way,  a check  was  made  on  the  assumption 
that  the  mixed  mean  water  temperature  remained  constant  through  the  test 
section.  This  tank  was  also  heavily  insulated. 

5.  Temperature  Control  Tank 

Immediately  downstream  of  the  mixing  chamber  was  the  temperature 
control  tank  which  was  used  to  remove  energy  from  the  working  fluid  as  it 
flowed  around  the  loop.  This  tank  was  open  and  cylindrical  in  shape  having 
a diameter  of  38.1  cm  (15  inches)  and  normally  holding  about  26  liters  of 
water.  The  tank  was  also  insulated  with  polyurethane  foam.  It  was  found 
that  the  energy  transferred  to  the  water  as  it  flowed  through  the  forced 
circulation  pump  exceeded  the  heat  transfer  from  the  water  as  it  flowed 
past  the  ice  layer.  Therefore,  a cooling  probe  was  placed  in  the  temper- 
ature control  tank  to  remove  this  excess  energy  and  thus  assure  that  the 
water  would  remain  at  essentially  constant  temperature  throughout  a given 
test.  The  cooling  probe  was  supplied  refrigerant  by  a mechanical 
refrigeration  unit  manufactured  by  FTS  Systems  Inc..  The  refrigeration 
system  could  be  operated  in  one  of  two  modes,  constant  cooling,  or 
thermostatically  controlled  cooling  in  which  temperature  sensing  was 


achieved  through  a copper-constantan  thermocouple  located  in  the  tank. 
Before  a run  began,  the  refrigeration  system  was  used  to  cool  the  entire 
flow  loop  down  to  the  desired  temperature  while  a small  auxiliary  pump 
provided  the  necessary  circulation.  A secondary  function  of  the  tempera- 
ture control  tank  was  that  of  filtering  the  water  as  it  flowed  through 
the  system.  Several  layers  of  fine  nylon  screen  covered  the  tank  exit 
(suction  side  of  the  forced  convection  pump)  and  thus  removed  any  foreign 
material  from  the  loop  which  may  have  been  introduced  through  the  open 
tank. 

6.  Forced  Circulation  Pump 

The  water  flow  in  the  loop  was  driven  by  the  main  pump  located 
somewhat  downstream  of  the  temperature  control  tank.  The  pump,  manu- 
factured by  the  Allis  Chalmers  Corp. , was  of  the  centrifugal,  closed 
impeller  design.  It  was  driven  by  a 3/4  horsepower  close-coupled 
electric  motor.  It  was  found  that  this  combination  was  capable  of 
providing  volumetric  flow  rates  of  up  to  1.76  liters/sec  (28  GPM)  and 
flow  velocities  of  1.2  m/sec  (3.9  ft/sec)  corresponding  to  a flow 
Reynolds  number  in  the  test  section  of  46,000.  The  flow  rate  through 
the  system  was  controlled  by  means  of  a throttling  valve  on  the  pump 
discharge. 


7.  Flow  Meter 

Flow  measurement  was  achieved  by  passing  the  flow  through  a 250  mm, 
high  capacity,  liquid  rotameter  manufactured  by  the  Ametek  Corporation. 

The  stainless  steel  float  was  non-viscosity  compensating  but  did  not 
alter  the  overall  accuracy  of  the  flow  measurements  because  the  water 
temperature  range  used  during  steady-state  and  transient  runs  represented 
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a fairly  small  variation  in  viscosity.  A discussion  of  rotameter 
calibration  is  included  in  the  following  section  on  temperature . flow, 
and  ice  thickness  measurements. 

C.  Tamper at ufe » Flow,  and  Ice  Thickness  Measurements 

As  mentioned  previously,  the  two  primary  temperature  measurements 

made  in  the  experimental  investigation  were  the  cooling  plate  surface 

temperature  and  bulk  water  temperature.  The  plate  surface  temperature 

in  the  immediate  vicinity  of  the  ice  thickness  measurement  was  used  in 

forming  the  dimensionless  parameter,  St  * c(T  -T  )/C  , while  the  bulk 

m w 

water  temperature  was  used  indirectly  in  determining  q" , the  heat  flux 

c 

at  the  ice-water  interface  required  to  compute  the  dimensionless  heat 

flux  parameter,  q*  » q"  S /fcp a . 

c o 

1 . Temporal uro  Measurement 

Cooling  plate  surface  temperatures  wore  obtained  from  two  thermo- 
couples imbedded  in  the  plate  Just  below  the  plate  surface,  immediately 
upstream  and  downstream  of  the  location  where  ice  thickness  measurements 
were  made.  The  copper -const ant  an  thermocouple  wires  were  0.072  mm 
(0.003  inch)  in  diameter  and  wore  insulated  by  a teflon  coating.  The 
copper-const antan  wire  combination  was  chosen  for  the  plate  thermocouples 
because  wires  made  from  these  materials  are  normally  of  a high  degree  of 
homogeneity  and,  therefore,  temperature  gradients  existing  along  t lie  wires 
do  not  cause  the  generation  of  extraneous  T.MF's  within  the  thermocouple 
circuit..  The  two  wires  of  each  thermocouple  were  laid  side-by-side  in 
two  shallow  grooves  which  were  machined  into  t lie  surface  of  the  copper 
cooling  plate  as  shown  in  Figure  6.  F.ach  groove  was  0.7b  mm  deep  and 
1.02  mm  in  width  and  extended  from  the  edge  of  the  plate  to  slightly  over 
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2.54 


Figure  6. 


Diagram  Illustrating  the  Method  of  Plato 
f>ur race  Temperature  Measurement 
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halfway  across  it,  perpendicular  to  the  direction  of  flow.  At  the  end 
of  each  groove  the  "floor"  of  the  groove  sloped  downward  to  meet  the 
plane  of  the  plate  surface  in  a gradual  manner,  see  Figure  6.  This  con- 
figuration allowed  the  junction  of  the  thermocouple  to  lie  very  near  to 
the  plate  surface.  The  junction  itself  was  made  by  a small  spot  weld  at 
the  wire  ends.  The  welded  junction  was  soldered  to  the  plate  at  the  end 
of  the  groove.  The  excess  solder  was  filed  off  level  with  the  plate 
surface  and  then  polished  smooth  with  a fine  abrasive.  The  remainder  of 
the  open  groove  between  the  thermocouple  junction  and  plate  edge  was 
filled  in  with  cement  and  sanded  smooth  and  level  with  the  plate  surface. 
The  problem  normally  associated  with  accurate  surface  temperature 
measurement,  that  of  maintaining  thermal  equilibrium  between  the  thermo- 
couple junction  and  the  surface,  was  minimized  in  this  instance  by  the 
high  thermal  conductivity  of  the  copper  from  which  the  cooling  plate  had 
been  fabricated.  Measured  heat  fluxes  through  the  plate  in  the  direction 
perpendicular  to  the  plate  surface  indicate  that  steady  state  temperature 
gradients  in  the  plate,  in  this  normal  direction,  were  a maximum  of  about 
1.5°  C/cm.  Therefore,  errors  in  surface  temperature  measurement  were  not 
considered  to  be  large  because  the  total  temperature  drop  through  the 
cooling  plate  was  less  than  1°  C.  However,  in  order  to  verify  that  the 
overall  error  involved  in  the  plate  surface  temperature  measurement  was 
indeed  small,  the  following  check  was  used.  A thin  ice  layer  was  grown 
over  the  entire  cooling  plate  surface  by  adjusting  the  convective  heat 
flux  at  the  ice-water  interface  through  manipulation  of  the  water  flow 
rate  in  the  channel.  The  convective  flux  was  then  increased  slowly  by 
increasing  the  water  flow  in  small  increments  until  the  temperature  of 
the  extreme  upstream  edge  of  the  plate  rose  above  the  fusion  temperature. 
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This  resulted  in  the  downstream  portion  of  the  plate  being  covered  by  a 
thin  patch  of  ice.  As  flow  rate  was  increased  still  further,  the  forward 
edge  of  the  ice  patch  moved  downstream  toward  the  surface  temperature 
measurement  positions.  As  the  edge  of  the  patch  reached  and  crossed 
over  these  locations,  the  thermocouple  junction  outputs  were  monitored 
on  a strip  chart  recorder  so  that  the  measured  surface  temperature  could 
be  compared  with  the  actual  surface  temperature  (assumed  to  be  0°C,  the 
fusion  temperature)  at  the  precise  instant  at  which  the  ice  patch 
"uncovered"  the  surface  thermocouples.  This  procedure  was  repeated 
several  times  for  each  of  the  two  surface  thermocouples.  The  results 
showed  that  each  thermocouple  measured  temperature  consistently  to 
within  0.25°C  of  0°C.  The  assumption  made  here  that  the  actual  surface 
temperature  was  0°C  is  supported  by  the  findings  of  Siegel  and  Savino  (7) 
in  which  the  temperature  of  the  ice-water  interface  under  similar 
experimental  conditions  was  always  found  to  be  within  O-l^C  of  0°C. 

Each  wire  of  the  two  copper-constantan  thermocouples  was  run  out 
of  the  flow  channel  through  grooves  in  the  cork  gasket  which  was 
installed  between  the  upper  surface  of  the  flow  channel  and  the  copper 
plate.  Both  thermocouples  were  referenced  at  0°C  with  a "Frigistor" 
automatic  ice  point  reference  junction.  The  two  plate  surface  thermo- 
couple outputs  could  be  read  by  strip  chart  recorder  or  through  an 
integrating  digital  voltmeter,  or  both. 

Bulk  water  temperatures  were  measured  with  chromcl-constantan 
thermocouple  probes  located  25.4  cm  (10  inches)  upstream  of  the  cooling 


plate  leading  edge  and  also  downstream  of  the  test  section  in  the  mixing 
chamber.  The  upstream  probe  was  of  the  exposed  junction  type  and  extended 
through  the  floor  of  the  flow  channel  so  that  the  junction  was  located 


approximately  at  the  geometric  center  of  the  channel  cross  sectional  flow 
area.  The  centerline  fluid  temperature  indicated  by  the  probe  was  taken 
to  be  identical  with  the  bulk  fluid  temperature  at  the  upstream  end  of 
the  test  section  because  of  the  heavy  insulation  which  covered  the  flow 
channel  and  inlet  piping.  Errors  in  bulk  water  temperature  measurement 
due  to  heat  conduction  along  the  stainless  steel  thermocouple  probe 
sheath  were  considered  negligible  for  several  reasons.  First  of  all,  the 
overall  temperature  difference  between  the  surrounding  atmosphere  and 
the  flowing  water  was  always  less  than  15°C;  thus,  there  was  little 
driving  potential  for  heat  conduction  along  the  probe  sheath.  Also,  the 
base  of  the  probe  was  threaded  into  the  plexiglass  flow  channel  which  in 
turn  was  heavily  insulated.  Therefore,  the  flow  channel  and,  hence,  the 
base  of  the  probe  could  not  have  been  at  a temperature  significantly 
different  from  that  of  the  water,  reducing  even  further  the  potential 
for  heat  flow.  Finally,  the  water  flow  velocities  through  the  flow 
channel  (0.35  m/sec.  to  0.60  m/sec.)  would  produce  heat  transfer 
coefficients  over  the  probe  of  such  a magnitude  that  heat  conducted  to 
the  probe  junction  (possibly  through  the  thermocouple  lead  wires)  would 
be  efficiently  transferred  to  the  flowing  water  through  a very  small 
temperature  difference. 

The  downstream  probe  was  used  primarily  to  verify  that  the  bulk 
water  temperature  remained  essentially  constant  as  it  passed  through  the 
test  section.  This  chromel-constantan  probe  was  of  the  grounded-junction 
type  and  extended  about  7.6  cm  (3  inches)  into  the  mixing  chamber.  For 
reasons  similar  to  those  outlined  above,  conduction  error  was  also 
considered  negligible  in  the  temperature  measurements  made  with  this 
probe.  Both  bulk  water  temperature  thermocouple  probes  were  referenced 
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at  0°C  with  the  Frigistor  automatic  ice  point  reference.  The  outputs  of 
each  probe  could  be  read  from  a digital  voltmeter  or  continuously 
monitored  by  strip  chart  recorder.  The  maximum  uncertainties  involved 
in  the  measurement  of  bulk  and  plate  surface  temperatures  (including 
uncertainties  due  to  the  limited  reading  precision  of  the  strip  chart 
recordings)  were  estimated  as  0.5°C  and  0.2°C  for  the  plate  and  bulk 
water  temperature  measurements,  respectively. 

2 . Flow  Measurement 

Water  flow  rate  measurement  in  the  test  loop  was  aclv'eved  with  a 
standard  liquid  rotameter  with  a 250  mm  meter  tube.  The  rotameter  float 
was  of  the  non-viscosity  compensating  or  streamlined  type;  however,  this 
did  not  significantly  affect  flow  measurement  accuracy  because  bulk  water 
temperatures  varied  only  over  a small  range  during  the  experiments  (less 
than  10°C).  The  rotameter  was  calibrated  directly  with  tap  water  at 
temperatures  within  the  range  used  in  both  steady-state  and  transient 
experiments.  The  36-point  calibration  yielded  a calibration  curve  which 
fit  the  36  data  points  with  a root  mean  square  deviation  of  0.43  liters/min. 
for  volumetric  flo  ; between  18.9  liters/min.  and  104  liters/min.  The 
accuracy  with  which  the  flow  rate  could  be  measured  depended  also  upon 
the  degree  of  precision  to  which  the  instrument  could  be  read.  This  was 
not  a limiting  factor  in  the  present  study  because  the  flute-guided  meter 
tube  allowed  the  float  position  to  be  estimated  to  within  + 1 mm  (+  0.40 
liters/min.)  under  most  circumstances.  Some  degree  of  float  oscillation 
did  occur  at  low  flow  rates;  under  these  conditions,  the  mean  flow  rate 
was  considered  to  correspond  to  the  midpoint  of  the  total  range  of  float 
oscillation.  The  absolute  maximum  uncertainty  in  the  flow  measurements 
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was  estimated  as  + 0.87  liters/min.  (corresponding  to  the  largest  devia- 
tion of  the  data  from  the  fitted  calibration  curve)  which  represents 
three  percent  of  the  lowest  flow  rate  encountered  in  any  part  of  the 
experiment . 

3.  loo  Layer  Thickness  Measurement 

Accurate  measurement  of  the  thickness  of  the  ice  layer  which  formed 
on  the  copper  cooling  plate,  both  at  steady-state  and  under  transient 
conditions,  was  essential  in  the  experimental  investigation.  There  are 
a number  of  ways  in  which  the  measurement  can  be  made,  all  of  which  fall 
into  basically  two  categories:  optical  methods,  and  mechanical  probe 
methods.  A short  discussion  of  the  measurement  techniques  of  each  type 
1 o 1 low. 


a.  Opt  tea)  Measurement  Teoliti  f<jues 
The  Ice  layer  thickness  can  he  measured  by  photographing  the  ice 
layer  in  a direction  normal  to  the  water  flow  and  in  t lie  plane  of  the 
ice-water  interface  as  si  town  in  Figure  7.  The  photographs  taken  in  this 
manner  are  then  enlarged  and  the  ice  layer  thickness  is  determined  from 
the  enlargements.  This  method  was  one  of  those  used  by  Siegel  and  Savino 
(7)  in  an  experimental  study  similar  to  tin*  present  one.  From  the  figure, 
it  is  seen  that  this  method  of  ice  thickness  measurement  is  useful  only 
when  the  ice  layer  is  ol  fairly  constant  thickness  across  the  channel. 

The  photographic  method  yields  measurements  ol  the  thickest  portion  ol 
the  lee  layer  and  cannot  he  used  to  measure  a "point  thickness"  at  some 
specific  location  on  the  plate.  Other  factors  which  limit  the  accuracy 
ol  the  photographic  technique,  according  to  Siegel  and  Savino,  are 
optical  distortion  efleets,  optical  illusions,  dltt  lenities  in  obtaining 
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exact  camera  alignment  with  the  plane  of  the  ice  interface,  and  limited 
image  magnification.  Optical  distortions  resulted  from  the  temperature 
gradient  in  the  thermal  boundary  layer  immediately  adjacent  to  the 
Interface,  and  optical  illusions  arose  due  to  the  transparency  of  the 
ice  layer  surface.  Image  magnification  was  a problem  despite  the  use  of 
a 70  mm  camera  in  the  Siegel  and  Savino  study.  Transient  ice  growth 
measurement  by  the  photographic  technique  involved  special  difficulties 
in  that  the  camera  alignment  with  the  ice-water  interface  had  to  be 
continuously  adjusted  throughout  the  transient. 

A second  optical  method  of  ice  thickness  measurement,  also  dis- 
cussed by  Siegel  and  Savino,  is  achieved  by  visual  means  rather  than 
photographic.  This  measurement  technique  is  also  illustrated  in  Figure  7, 
in  which  a small  telescope  mounted  on  a cathetometer  is  used  to  visually 
locate  the  ice-water  interface.  With  this  method,  grids  of  horizontal 


lines  spaced  at  a known  distance  apart  are  placed  on  both  walls  of  the 
flow  channel.  When  the  interface  is  sighted  will)  the  telescope  in  the 
same  plane  as  two  corresponding  grid  lines  on  either  side  of  the  channel, 
t ho  ice  layer  thickness  is  simply  read  from  the  grids.  According  to 
Siegel  and  Savino,  this  technique  was  found  to  suffer  from  many  of  the 
same  disadvantages  as  the  photographic  method;  i.e.  , optical  distortion 


due  to  temperature  gradients  in  the  boundary  layer  and  difficulties  in 
locating  the  transparent  ice  surface. 

Neither  of  the  opt  leal  methods  described  above  are  ideal  for 
measurement  of  solid  layer  thicknesses  at  a specific  point  on  the  cooling 
plate.  The  mechanical  and  temperature  probe  methods  described  below  do 
not  have  this  limitation. 
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b.  Mechanical  and  Temperature  Probe  Techniques 
The  mechanical  probe  method  of  ice  thickness  measurement  consists 
of  introducing  a small  probe  into  the  test  section  through  the  channel 
wall  opposite  the  ice  layer  and  raising  the  probe  until  it  contacts  the 
layer  surface,  thus  yielding  a direct  measurement  of  the  ice  layer  thick- 
ness at  a given  point  on  the  cooling  plate.  A variation  of  this  procedure 
was  used  by  Siegel  and  Savino  in  which  the  probe  was  fabricated  in  the 
form  of  a fine  gauge,  forked  thermocouple  probe,  see  Figure  8.  The 
thermocouple  wire  extending  between  the  fork-like  supports  was  0.076  mm 
(0.003  inch)  in  diameter  with  a butt-welded  junction  at  the  center.  The 
probe  shaft  was  mounted  in  a micrometer  head  so  that  it  could  be  raised 
or  lowered  in  fine  increments.  In  measuring  the  ice  layer  thickness,  the 
probe  was  made  to  approach  the  ice  surface  slowly  while  monitoring  the 
output  of  the  probe  thermocouple  on  a strip  chart  recorder.  When  the 
thermocouple  output  corresponded  to  0°C,  the  fusion  temperature,  the  probe 
was  assumed  to  be  in  contact  with  the  ice  layer  surface  and  the  thickness 
was  then  read  from  the  micrometer  dial.  According  to  Siegel  and  Savino, 
this  temperature  probe  method  of  ice  thickness  measurement  proved  to  be 
most  accurate  and  was  reported  to  yield  thickness  measurements  correct  to 
within  0.025  mm  (0.001  inch).  In  addition  to  the  temperature  probe,  a 
purely  mechanical  probe  was  also  employed  in  the  Siegel  and  Savino 
experiments  to  measure  ice  layer  thickness  under  steady-state  conditions. 
This  probe  consisted  of  a thin  metal  rod  which  was  mounted  in  a micrometer 
head  in  a manner  similar  to  the  temperature  probe.  The  rod  was  moved 
toward  the  interface  until  it  was  visually  determined  that  the  probe  had 
contacted  the  ice  layer,  at  which  time  the  ice  thickness  was  read  from  the 
micrometer.  It  was  not  stated  by  the  authors  to  what  degree  of  accuracy 
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Figure  8.  Illustration  of  the  Thermocouple  Probe 
Ice  Thickness  Measurement  Technique 
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the  lee  thickness  could  he  measured  by  this  method . Both  the  temperature 
and  mechanical  probe  methods  described  here  raise  significant  questions' 
will*  respect  to  their  accuracy  In  the  experimental  determination  of  ice 
layer  thickness.  Because  a similar  type  of  ice  probe  was  used  in  the 
present  experimental  investigation,  the  advantages,  disadvantages,  and 
uncertainties  Inherent  in  their  use  will  be  examined  in  some  detail. 

c . Unct-rt -itut v in  Me ch ; i nlcal  and  Temperature  Probe  Techniques 
and  Sirepl ified  Error  Analysis 

in  order  to  accurately  measure  ice  layer  thickness  by  either 
mechanical  or  temperature  probe  methods , the  presence  of  the  probe  must 
have  a minimal  effect  on  the  ice  layer  surface  as  the  probe  is  moved  into 
the  vicinity  of  the  layer.  Otherwise,  the  measured  thickness  will  not 
represent  the  true,  undisturbed  ice  layer  thickness  for  the  prevailing 
conditions  (plate  surface  temperature  and  convective  heat  flux  at  the 
interface).  The  phenomenon  under  consideration  here  is  somewhat  analogous 
to  the  so-called  "loading  error"  which  may  occur  in  electrical,  transducer- 
typo  measurement  devices.  The  fundamental  principal  which  governs  in  both 
instances  is  that  the  measurement  process  will  inevitably  alter  the 
quantity  which  is  to  be  measured,  see  Reference  17.  Errors  of  this  type 
arise  in  mechanical  probe  thickness  measurements  primarily  because  the 
probe  must  be  of  finite  dimension  and  therefore,  when  brought  into  the 
immediate  vicinity  of  the  ice-water  interface,  will  have  some  effect  on 
the  convective  heat  transfer  phenomena  occurring  there.  The  effect 
probably  can  be  attributed  to  a local  flow  distortion  caused  by  the  probe's 
presence  which  in  turn  perturbs  the  local  convective  heat  flux  at  the 
ice-water  interface  from  Its  "normal"  (in  the  absence  of  the  probe)  value. 
This  flow  distortion  effect  will  always  be  present  to  some  degree  and  its 
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magnitude  will  probably  depend  heavily  on  the  probe  size  and  shape  in 
addition  to  the  flow  conditions  at  the  interface.  In  the  present  experi- 
mental study,  the  ice  thickness  always  tended  to  decrease  very  slightly 
in  the  extended  presence  of  the  probe,  thus  implying  a slightly 
increased  heat  transfer  coefficient  at  the  interface.  To  minimize  the 
distortion  of  the  ice  layer  due  to  the  probe,  steps  may  be  taken  in  the 
design  of  the  probe  and  also  in  the  experimental  technique. 

Ice  thickness  measurement  by  the  temperature  probe  method  similar 
to  that  used  by  Siegel  and  Savino  is  subject  to  an  even  wider  range  of 
experimental  errors.  In  addition  to  the  requirement  that  the  probe  should 
not  disturb  the  ice-water  interface,  the  basic  principal  of  the  measure- 
ment technique  requires  that  the  probe  tip  or  thermocouple  junction  be 
very  nearly  in  thermal  equilibrium  with  its  surroundings  when  the 
junction  is  immediately  adjacent  to  the  ice  layer  surface.  In  other 
words,  the  conduction  error  must  be  kept  small  if  accurate  ice  thickness 
measurements  are  to  be  made.  In  order  to  gain  a better  qualitative 
understanding  of  the  immediate  conduction  error  problem  and  to  identify 
the  effects  of  probe  size,  bulk  water  temperature,  flow  conditions,  etc., 
a simplified  analysis  was  carried  out  for  a cylindrical  probe  projecting 
into  the  thermal  boundary  layer  adjacent  to  the  solid-liquid  interface. 

Consider  the  situation  shown  in  Figure  9.  The  cylindrical  probe 
of  diameter  "d"  projects  perpendicularly  into  the  thermal  boundary  layer 
and  virtually  contacts  the  ice-water  interface.  For  the  purposes  of  this 
analysis,  it  was  assumed  that  heat  is  transferred  to  or  from  the  probe 
tip  by  forced  convection.  The  water  which  flows  past  the  probe  and  over 
the  interface  is  assumed  to  flow  under  the  influence  of  a negligibly  small 
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pressure  gradient  and  with  thermal  and  hydrodynamic  boundary  layer 
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thicknesses  of  comparable  dimension  (<5^  = 5 *“  5)  • Further  simplifying 

assumptions  made  in  the  analysis  are  listed  below. 

1.  Conduction  heat  transfer  within  the  probe  is  one-dimensional 
and  in  the  direction  perpendicular  to  the  ice  layer  surface 
(the  probe  has  an  effective  thermal  conductivity  in  this 
direction  of  "X"). 

2.  At  the  outer  edge  of  the  thermal  boundary  layer  Z = <5 , the 

probe  temperature  corresponds  to  that  of  the  freestream,  T . 

00 

3.  The  fluid  temperature  in  the  boundary  layer  (T  ) varies  with 

r 

Z from  the  melting  temperature  (T  ) at  Z = 0 to  the  free- 

m 

stream  temperature  (T^)  at  Z - 6 . The  specific  dependence 
of  Tp  on  Z will  be  discussed  after  developing  the  model. 

4.  In  the  interval  0 < Z £ 6 , a uniform,  average  film  coefficient 
exists  over  the  probe  surface  of  magnitude  h . 

5.  The  tip  of  the  probe  transfers  heat  to  or  from  the  fluid  via 
the  effective  film  coefficient  h^  which  is  based  upon  the 

temperature  difference  (T  _ - T ). 

Z=u  m 

The  validity  of  the  uniform  film  coefficient  assumption  probably 
depends  upon  the  magnitude  of  the  actual  variation  of  the  conductance  in 
the  interval  0 < Z 6 . If  this  variation  is  small,  the  assumption  of 
uniformity  in  the  film  coefficient  should  not  introduce  significant  error 
Into  the  analysis.  It  is  shown  below  that  the  local  film  coefficients  at 
any  two  points  on  a small  diameter  probe  in  a low  speed  flow  are  of  the 
same  order  of  magnitude. 

Consider  the  crossflow  of  a liquid  over  the  cylindrical  probe  shown 
in  Figure  10.  For  10  < Red  = Vd/v  < 1000,  the  boundary  layer  over  much 


of  the  probe  periphery  is  laminar  with  small  vortices  being  shed  in  the 

vicinity  of  the  downstream  stagnation  point  (Reference  18).  Within  this 

range  of  Red,  the  character  of  the  boundary  layer  flow  and  therefore  the 

heat  transfer  phenomena  on  the  cylinder  do  not  change  significantly  with 

Red.  Thus,  it  is  expected  that  Nud  for  this  Red  range  is  not  a strong 

function  of  Red.  For  the  flow  conditions  encountered  in  the  experimental 

apparatus  described  earlier,  the  maximum  value  of  Red  was  always  less 

than  200  based  on  the  freestream  velocity.  From  the  above  considerations, 

it  can  be  concluded  that  a largely  laminar  boundary  layer  exists  over  the 

probe  for  all  values  of  Z between  Z = 0 and  Z = 6 providing  an 

essentially  uniform  film  coefficient  over  the  entire  probe.  The  variation 

of  viscosity  between  Z = 0 and  Z = 6,  also  reduces  the  variation  in  h 

along  the  probe  shaft.  As  Prandtl  number  increases  at  constant  Red,  Nud 

1/  3 

increases  at  a rate  of  approximately  Pr  . For  example,  consider  the 
correlation  due  to  McAdams  (18)  for  liquids  in  crossflow  over  cylinders, 

Nud  = 1. l(Red)n  Pr0,31  . (22) 

It  was  determined  from  this  equation  and  from  typical  property  values  and 
flow  conditions  encountered  in  the  present  experiment  that  Nud  varied  by 
a factor  of  about  2 between  two  positions  on  the  probe  for  which  Red 
varied  by  an  order  of  magnitude.  Therefore,  the  uniform  film  coefficient 
assumption  is  reasonable,  at  least  for  the  level  of  analysis  to  be 
carried  out  here. 

With  these  assumptions,  an  energy  balance  on  a differential 
element  of  the  probe  yields  the  following  simple  system: 
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Ad“  d T - „ . 

~ 4 J ~ hd(T-Tp)  = 0 0 < Z < 6 

dZ~ 


A f - ht<T-T«>  ■ 0 Z - 0 . 


By  introduction  of  the  dimensionless  quantities. 


A = 


T-T 


T -T 


and 


(23a) 


(23b) 


(23c) 


and  substitution  into  Equations  (23),  the  following  dimensionless  system 
results: 

A"  - y2A  = -Y2Ap  0 < Y < 1 (24a) 


A = 1 


Y = 1 


(24b) 


1 


A’  - 3A  * 0 


Y = 0 . 


(24c) 


The  parameters  y and  g are  given  by, 


2 - 4h  x2 


(25) 


V 


A * 


(26) 


while  Ap  is  the  dimensionless  local  fluid  temperature  in  the  boundary 
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layer. 


T -T 
1 m 

T -T 
05  m 


(27) 


2 

For  a given  boundary  layer  thickness,  Y and  8 indicate  the  magnitudes 

of  the  film  coefficients  over  the  probe  shaft  and  probe  tip,  respectively. 

The  deviation  from  zero  of  A(Q)  , the  dimensionless  probe  tip  temperature, 

is  a direct  indication  of  the  amount  by  which  the  tip  temperature  differs 

from  the  melting  temperature,  T , and  hence,  is  a measure  of  the  conduc- 

m 

tion  error  which  would  affect  the  accuracy  of  the  measurement.  The  general 
solution  to  the  system  (24)  is  given  by  Equation  (28). 


A(Y) 


Asinh  Y V + Bcosh  Y Y 


- Y 


jApSinh  Y(Y-n) 
o 


dn 


(28) 


where  the  constants  A and  B are  obtained  from  conditions  (24b)  and 
(24c): 


A = [f-  + 0 [ A sinh  Y(l-rl)  dn]/[~;  sinh  Y + cosh  y]  (29) 

i ) t Y 

. o 

and 

B = [1  + Y | ApSinh  Y ( 1— n)  dq]/[^  sinh  Y + cosh  Y]  • (30) 

o 

The  temperature  profile  in  the  boundary  layer  T (2)  has  a 

r 

substantial  effect  on  the  probe  tip  temperature.  Therefore,  three 
temperature  profiles  will  be  considered  here,  two  of  which  correspond  to 
the  commonly  encountered  hydrodynamic  regimes  of  laminar  and  turbulent 
boundary  layer  flow.  The  temperature  variations  in  the  boundary  layer 
were  taken  to  be  the  same  as  those  commonly  employed  in  the  thermal 


a 


analysis  of  two-dimensional  boundary  layers  via  the  energy  integral 
equation.  Each  profile  can  be  expressed  by: 
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af  = 


C1Y 


C V 

V 


(31) 


where  the  values  of  the  constants  C^,  C and  k in  Equation  (31)  for 
the  linear,  laminar,  and  turbulent  temperature  profiles  are  shown  in 
Table  1 (Reference  19  and  Reference  14).  The  profiles  are  also  plotted 
in  Figure  11. 

Substitution  of  the  linear  profile  into  Equations  (28),  (29),  and 
(30)  results  in  the  following  simple  relationship  between  A and  Y : 


A(Y)  = Y + 


sinhy ( 1-Y) 
ycoshy  + gsinhy 


(32) 


The  dimensionless  probe  tip  temperature  becomes: 


A(0) 


1 

y co thy  + 3 


(33) 


Note  that  A(0)  tends  toward  zero  when  y and/or  3 become  large.  Thus 
it  is  possible,  at  least  in  the  case  of  a linear  temperature  profile,  to 
decrease  the  difference  between  the  probe  tip  temperature  and  the  melting 
temperature  by  increasing  the  film  coefficients  over  the  probe  shaft 
and/or  probe  tip. 

For  the  laminar  and  turbulent  temperature  profile,  A(Y)  was 
obtained  through  numerical  integration  of  Equation  (28)  via  Simpson's 
3/8  rule.  The  results  are  shown  in  Figures  12  and  13. 
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Figure  11.  Linear,  Laminar,  and  Turbulent  Boundary 
Layer  Temperature  Profiles 


Figure  12.  Dimensionless  Probe  Temperature  A 
Versus  Y 
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Figure  12  is  a plot  of  the  dimensionless  probe  temperature.  A, 
as  a function  of  Y for  all  three  temperature  profiles  with  6 as  a 
parameter  and  y = 10  . From  this  figure,  there  is  seen  to  be  a sub- 
stantial variation  in  probe  temperature  due  only  to  differences  in  the 
boundary  layer  temperature  profile,  all  other  conditions  held  constant. 
Also,  the  effect  of  increased  film  conductance  at  the  probe  tip  or 
increased  thermal  boundary  layer  thickness  (increased  (?)  serves  to 
decrease  the  dimensionless  probe  temperature  for  small  values  of  Y . 

For  Y values  approaching  unity,  the  convective  bound. i-"  condition  at 

Y = 0 has  little  effect  and  A is  seen  to  approach  A_,  the  local 

F 

boundary  layer  fluid  temperature  (compare  with  Figur  ?.  This 
behavior  for  die  laminar  and  turbulent  boundary  layeL  Vs  parallels 

that  found  in  Equation  (32)  for  the  linear  profile. 

The  critical  quantity  with  regard  to  thermocouple  probe  measure- 
ments of  ice  layer  thickness  is  A(0),  the  dimensionless  probe  tip 
temperature.  A zero  value  of  A(0)  indicates  that  the  probe  tip  is  at 
the  same  temperature  as  the  ice-water  interface.  In  reality,  however, 
A(0)  can  only  approach  zero  if  the  probe  tip  remains  in  the  position 
shown  in  Figure  9.  The  tip  temperature  can  be  made  equal  to  the  melting 
temperature  only  by  advancing  the  warm  probe  tip  into  the  ice  layer,  thus 
distorting  the  layer  surface  and  yielding  false  thickness  readings.  The 
magnitude  by  which  A(0)  deviates  from  zero  when  the  tip  is  immediately 
adjacent  to  the  undeformed  layer  is  shown  in  Figure  13.  The  significance 
of  the  boundary  layer  temperature  profile  is  readily  obvious  from  the 
curves  in  this  figure.  For  all  values  of  the  parameter  £?  shown,  the 
values  of  A(0)  indicated  Cor  the  laminar  and  turbulent  boundary  layer 
temperature  profiles  differ  by  nearly  an  order  of  magnitude  for  values  of 
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the  parameter  Y greater  than  about  20.  For  low  values  of  y,  the 
curves  for  the  laminar  and  turbulent  profiles  converge  and  eventually 
meet  at  Y » 0 . The  coincidence  of  the  curves  at  y ■ 0 derives  from 
the  fact  that  the  probe  is  perfectly  insulated  for  y » 0 and  A(0)  is 
therefore  independent  of  temperature  variations  in  the  boundary  layer. 

A somewhat  surprising  feature  of  the  curves  in  Figure  13  is  that  A(0) 
is  not  necessarily  a monotonically  decreasing  function  of  y . For  large 
values  of  the  probe  tip  film  coefficient  h (large  values  of  (3) , A(0) 

increases  for  a time  with  increased  h (increased  y)  and  then  decreases 
toward  zero  for  larger  values.  This  effect  is  seen  to  be  especially 
pronounced  in  the  case  of  a turbulent  boundary  layer  temperature  profile. 

In  fact,  for  turbulent  profiles  and  values  of  P > 1,  conduction  error 
could  probably  be  more  effectively  reduced  by  insulation  of  the  probe 
shaft.  If  the  shaft  is  not  insulated,  then  to  decrease  the  conduction 
error,  y (or  h)  must  be  increased  to  a value  to  the  right  of  the 
maximum  in  the  A(0)  curve.  In  order  to  attach  some  physical  significance 
to  this  effect,  consider  Figure  14  In  which  the  probe  is  shown  schema- 
tically along  with  the  turbulent  fluid  temperature  profile.  For  y ■ 0 
(probe  shaft  insulated),  the  probe  is  in  thermal  communication  with  the 
surrounding  fluid  only  at  the  tip  (V  » 0)  via  the  dimensionless 
conductance  p and  the  temperature  profile  in  the  probe  is  linear.  For 
very  large  p,  the  tip  temperature  will  approach  A(0)  ■*  0 . If  y is 
increased  to  some  small,  non-zero  value,  the  tip  area  of  the  probe  (Y  » 0) 

will  bo  exposed  to  the  fluid  temperature  A.,  which  increases  very 

r 

rapidly  for  small  Y and  heat  will  be  transferred  to  the  probe  from  the 
fluid  by  convection  and  then  flow  toward  the  tip,  thus , increasing  the 
tip  temperature.  As  y is  increased  further  (moving  toward  the  right  on 
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Figure  14.  Sketch  of  Thermocouple  Probe  and 
Turbulent  Temperature  Profile 
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Figure  13),  this  effect  must  eventually  reverse  itself  because,  for 
infinite  shaft  conductance  (y  -*■<»),  the  temperature  profile  in  the 
probe  will  coincide  with  the  temperature  profile  in  the  fluid,  and  then 
A(0)  - 0 . 

It  must  be  concluded  from  the  preceding  analysis  that  ice  layer 
thickness  measurement  by  the  thermocouple  probe  technique  poses  many 
difficulties.  For  a laminar  boundary  layer  at  the  ice-water  interface, 
the  problems  could  probably  be  overcome  with  careful  probe  design  and 
use.  In  turbulent  flow,  the  values  of  the  film  conductance  parameters 
y and  f?  will  most  likely  be  of  the  same  order  as  in  laminar  flow. 

This  is  because  the  boundary  layer  over  the  shaft  of  a small  diameter 
probe  in  a low  speed  flow  will  be  mainly  laminar  regardless  of  surrounding 
flow  conditions.  If  the  laminar  and  turbulent  boundary  layer  thicknesses 
are  also  assumed  equal,  it  is  seen  from  Figure  13  that  conduction  error 
is  much  more  of  a problem  in  the  turbulent  case. 

d . The  Method  Used  in  this  Study 

It  has  been  shown  that  the  difficulties  in  obtaining  ice  layer 
thickness  measurements  with  the  thermocouple  probe  technique  are  sub- 
stantial. In  addition,  the  problems  wore  show  to  be  far  more  serious 
when  the  boundary  layer  flow  at  the  ice-water  interface  is  turbulent 
rather  than  laminar.  For  these  reasons,  a thermocouple  probe  technique 
was  not  used  in  this  study  to  measure  ice  layer  thickness.  All  steady- 
state  and  transient  thickness  measurements  were  made  by  means  of  a 
mechanical  probe  method  in  which  the  probe  tip  was  allowed  to  just 
contact  the  ice  layer.  Contact  between  probe  tip  and  ice  layer  was 
verified  visually  by  means  of  a largo  magnifying  lens  mounted  outside  of 
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the  test  section.  The  probe  used  was  a thermocouple  probe;  however,  the 
thermocouple  feature  of  the  probe  was  not  used  to  directly  measure  ice 
layer  thickness  but  only  to  aid  in  the  recording  of  transient  ice  thick- 
nesses. 

The  thermocouple  probe  was  stainless  steel  sheathed,  cylindrical, 
and  had  an  overall  diameter  of  0.25  mm  (0.01  inch).  The  thermocouple 
junction  was  made  from  0.04  mm  (0.0015  inch)  diameter  chromel  and 
constantan  wires  and  was  grounded  to  the  probe  sheath  at  the  tip.  The 
fine  gauge  probe  was  supported  by  two  larger,  concentric  stainless  steel 
tubes  which  projected  through  the  floor  of  the  flow  channel  below  the 
cooling  plate.  The  larger  tube  was  0.635  cm  (0.25  inch)  in  diameter  and 
15.24  cm  (6  inches)  in  length.  This  shaft  was  secured  to  the  slider 
block  of  a "Unislide"  micrometer  mechanism  mounted  underneath  the  test 
section.  The  "Unislide"  was  driven  by  a Starret  micrometer  head  with  a 
2.54  cm  (1  inch)  drive  range  accurate  to  +0.025  mm  (+0.001  inch).  The 
probe  is  shown  schematically  in  Figure  15  and  the  entire  assembly  in 
Figure  16.  In  order  to  minimize  local  ice  layer  distortion  in  the 
vicinity  of  the  probe,  the  flexible  probe  shaft  was  bent  almost  parallel 
to  the  ice  layer  surface  and  directed  into  the  water  flow  as  shown  in 
Figure  15.  In  this  configuration  a greater  length  of  t ho  probe  shaft 
was  immersed  in  the  thermal  boundary  layer  at  the  ice-water  interface 
and  thus,  the  conduction  of  heat  to  the  probe  tip  and  subsequently  to 
the  ice  surface  was  minimized.  The  ice  surface  was  distorted  a negligible 
amount  if  the  probe  was  advanced  toward  the  layer  quickly  and  immediately 
withdrawn  after  a measurement  was  made.  The  prolonged  presence  of  the 
probe  at  the  interface  did  cause  the  layer  to  experience  localized  melting. 
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Figure  15.  Mechanical-Thermocouple  Probe 


Tho  uncertainty  in  ice  layer  thickness  measurement  by  the  technique 
described  above  was  estimated  to  be  no  more  than  +0.05  ram  for  steady- 
state  thickness  measurements.  This  figure  is  based  upon  the  variation 
encountered  in  successive  thickness  measurements  made  while  the  cooling 
plate  temperature  and  bulk  water  temperature  remained  completely  static. 

D.  Steady-State  Heat  Flux  Measurements 

In  order  to  supply  the  finite  difference  model  with  accurate  values 
of  the  parameter  q*  for  comparison  with  the  experimental  results,  it 
was  essential  that  q'^,  the  interfacial  heat  flux,  be  known  for  a given 
flow  condition  and  bulk  fluid  temperature  in  the  channel.  Because  heat 
transfer  conditions  in  the  channel  did  not  correspond  to  any  simple 
solution  presently  existing  in  the  heat  transfer  literature,  it  was 
decided  that  the  most  reliable  and  direct  way  to  determine  the  convective 
film  coefficient  at  the  ice-water  interface  was  by  experiment. 

1.  Procedure 

In  the  experimental  procedure  used,  a steady-state  ice  layer  was 
formed  on  the  copper  cooling  plate  while  volumetric  water  flow  rate  and 
bulk  temperature  were  maintained  constant.  The  ice  layer  which  formed 
was  always  of  the  basic  shape  shown  in  Figure  17.  The  surface  of  the 
layer  was  smooth  and  level  in  the  direction  normal  to  the  flow  indicating 
little  or  no  cooling  plate  surface  temperature  variation  in  this  direction. 
Even  though  care  was  taken  to  build  up  the  ice  layer  slowly,  some  cracks 
would  usually  develop  in  the  layer  in  a direction  normal  to  the  plane  of 
the  cooling  plate  surface.  Because  the  cracks  always  formed  in  a direc- 
tion parallel  to  the  primary  heat  flow  path,  they  were  considered  to  have 
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a negligible  effect  on  the  heat  conduction  process  in  the  layer.  Other- 
wise, the  ice  was  clear  and  homogeneous  in  appearance. 

The  steady-state  ice  layer  thickness,  S,  was  used  to  calculate 

q",  the  interfacial  heat  flux,  from  the  relation: 
c 


(T  -T  ) 
m w 


(34) 


where  k is  the  ice  thermal  conductivity.  An  interfacial  film  coefficient, 
h,  was  defined  as 


h = 


(T, -T  ) 
b ra 


(35) 


where  T,  is  the  bulk  water  temperature.  Substituting  Equation  (35)  in 
b 

Equation  (34)  results  in 


(T  -T  ) 
m w 

(t, -t  ) 

b m 


(36) 


The  plate  temperature,  T , in  Equation  (36)  was  determined  from  the 

w 

outputs  of  the  two  plate  thermocouples  located  immediately  upstream  and 

dowstream  of  the  ice  thickness  measurement  position.  During  each  steady- 

state  run,  these  two  thermocouple  outputs  were  recorded  on  a dual-channel 

strip  chart  recorder.  T , in  Equation  (36) , was  given  by  the  average  of 

w 

these  two  temperatures.  Due  to  the  rather  close  proximity  of  these 
thermocouples  to  the  thickness  measurement  position  and  due  to  the  very 
high  thermal  conductivity  of  the  copper  plate,  this  linear  interpolation 
probably  did  not  introduce  significant  error  into  the  measurement  of  plate 
surface  temperature.  Also,  from  the  two  recorded  plate  temperatures,  a 
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plate  surface  tempo rat  tiro  gradient  was  calculated  as  an  experimental 
cheek  on  the  one-dimensional  heat  flow  assumption.  In  all  runs,  the 
temperature  gradient  along  the  plate  surface  In  the  direction  of  1 low 
was  less  than  1°  C/cm  and  on  the  average  was  less  than  half  this  magni- 
tude. The  temperature  gradients  through  the  ice  layer  (normal  to  the 
plate  surface),  on  the  other  hand,  were  always  in  the  range  150-!100°  C/em. 
From  these  results  and  the  fact  that  the  ice  layer  thickness  was  normally 
between  1 mm  and  4 mm,  it  would  seem  reasonable  to  assume  that  heat 
conduction  in  the  layer  was  very  nearly  one-dimensional.  T , the  fluid 
bulk  temperature,  was  indicated  by  a thermocouple  probe  upstream  of  t lie 
test  section.  This  temperature  was  also  recorded  on  a strip  chart 
recorder.  Several  ice  thickness  measurements  were  made  during  a single 
run  by  rapidly  moving  the  probe  tip  toward  the  interface  while  observing 
t lie  probe  tip  position  through  a magnifying  lens.  When  the  tip  contacted 
the  layer  surface,  the  micrometer  reading  was  noted  and  the  probe  quickly 
withdrawn  so  that  deformation  of  the  ice  layer  would  not  result.  Alter 
a run  was  completed,  the  plate  temperatures  and  the  results  ot  the 
thickness  measurements  were  examined  to  ascertain  which  measurements  were 
made  wh ( 1 e truly  steady-state  conditions  prevailed.  The  steady-state 
condition  was  indicated  by  plate  temperatures  which  changed  negligibly 
with  time  and  successive  thickness  measurement s which  showed  essentially 
iu>  so l i d 1 f l cat i on  or  melt  ing  ot  the  layer.  On  the  average,  a single  run 
c old  product'  1 to  5 thickness  measurements  taken  In  the  presence  of 

late  conditions.  Only  these  data  were  used  in  the  calculation  ot 


■i  it ui,  ditterencc  (T  -T  ),  an  Indication  of  the  sub- 
til w 


nil  , tanged  item  -2  1JC  to  -58*0  in  the  si  oadv-stat 
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experiments.  Subcooling  of  this  magnitude  indicates  that  the  thermal 
conductivity  variation  of  the  ice  could  be  appreciable.  Therefore,  the 
conductivity,  k,  in  Equation  (36)  was  computed  from  the  relation  given 
by  Ratcliff  (20): 


_)  = —UL 

V T ( K ) 


0.00615 


This  equation  is  valid  for  clear  ice  down  to  temperatures  of  -150°C.  The 

value  of  lc  in  Equation  (36)  was  evaluated  from  Equation  (37)  at  the 

temperature  (T  -T  )/2  which  represents  an  approximate  mean  temperature 
in  w 

in  the  layer. 

2 . Uncertainty  in  the  Prediction  of  Interface  Film  Coefficient 

In  performing  the  steady-state  measurements,  bulk  fluid  temperature 

and  water  flow  rate  were  chosen  so  that  the  resultant  steady-state  ice 

thickness  was  not  excessively  small.  Also,  care  was  taken  to  keep  the 

dry  ice-methanol  bath  in  a condition  which  would  assure  the  lowest 

possible  plate  temperatures.  These  precautions  were  taken  in  order  to 

minimize  error  in  the  measured  value  of  h . From  Equation  (36) , it  is 

evident  that  small  experimental  errors  in  plate  and  bulk  water  temperature 

measurement  and  ice  thickness  measurement  could  cause  significant  error 

in  the  computed  value  of  the  film  coefficient,  h,  if  the  magnitude  of 

the  errors  in  the  temperature  and  thickness  measurements  are  comparable 

to  the  quantities  (T  -T  ),  (T, -T  ),  and  S . The  maximum  possible 

m w b m 

error  est invited  for  h,  based  upon  the  uncertainties  in  the  measurements 

of  S,  T , and  T (discussed  in  the  previous  section)  did  not  exceed 
w b 

+7Z  for  most  of  the  steady-state  runs.  For  a few  runs  at  very  high 


interfacial  heat  flux  where  the  ice  layer  was  rather  thin  and  the  plate 

subcooling  (T  -T  ) small,  the  maximum  possible  error  in  h was 
m w 

estimated  to  be  as  much  as  +107,.  It  should  be  emphasized  that  these 
estimated  percentage  errors  in  h represent  upper  bounds  to  the  error 
and  most  probably  the  actual  error  in  the  prediction  of  film  coefficient 


was  much  smaller. 


E.  Transient  Ice  Growth  Measurements 


The  objective  in  the  transient  experiments  was  to  compare  experi- 
mental ice  layer  thickness  measurements  with  the  predictions  of  the  finit 
difference  model  developed  in  Chapter  II.  It  was  therefore  necessary  to 

obtain  T , q",  and  S as  functions  of  time  in  each  transient  run.  hue 
w c 

to  the  nature  of  the  mechanical  probe  technique  for  measuring  transient 
ice  layer  thickness,  only  transients  which  involved  growth  of  the  ice 
layer  were  studied.  Each  run  was  classified  as  to  the  method  by  which 
the  layer  growth  was  initiated.  There  were  basically  two  types  of 
transient  runs  made;  each  is  discussed  in  detail  below. 


1 • Plate  Temperature  Tra ns  ion ts 

In  this  first  type  of  transient,  a steady-state  ice  layer  was 
formed  on  the  cooling  plate  by  holding  plate  surface  temperature,  water 
flow  rate,  and  hulk  temperature  constant.  The  plate  surface  temperature 
in  this  steady-state  condition  was  maintained  above  -20°C  by  using  a 
reduced  amount  of  dry  ice  in  the  cooling  bin.  At  time  t = 0,  ice  layer 
growth  was  initiated  by  reducing  the  plate  temperature  while  water  flow 
rate  and  bulk  water  temperature  and  thus  interface  heat  flux  remained 
essentially  constant. 
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The  initial  thickness,  S^,  of  the  steady-state  ice  layer  was  not 
measured  directly  by  the  mechanical  probe  technique  but  instead  was 
computed  from  the  steady-state  plate  temperature  and  interfacial  heat 
flux,  the  interface  heat  flux  being  determined  from  the  bulk  water  temper- 
ature, the  water  flow  rate,  and  the  results  of  the  steady-state  heat  flux 
experiments.  This  procedure,  was  adopted  for  two  reasons.  First,  the 
initial  ice  layer  thickness,  S^,  for  the  plate  transient  runs  was  rather 
small  (less  than  1 mm)  so  that  a direct  thickness  measurement  with  the 

mechanical  probe  would  introduce  a +5%  to  +10 % error  into  the  S measure- 

— — o 

ment  (based  upon  +0.05  mm  uncertainty  in  the  thickness  measurement). 

Also,  due  to  the  difficulties  encountered  in  maintaining  a steady  plate 

temperature  with  a limited  amount  of  dry  ice  in  the  coolant  bin,  it  was 

very  difficult  to  obtain  an  accurate  measurement  of  with  the  probe. 

Under  these  circumstances,  prediction  of  the  initial  steady-state  ice 

thickness  S from  the  known  quantities  q"  and  T was  considered  to 
o c w 

be  the  most  reliable  method. 

Once  steady-state  had  been  achieved,  a large  quantity  of  dry  ice 
was  introduced  into  the  coolant  bin,  rapidly  lowering  the  plate  surface 
temperature  and  initiating  ice  growth.  The  growth  transient  normally 
lasted  approximately  A to  6 minutes  with  the  rate  of  growth  rapidly 
decreasing  with  time.  During  the  transient,  the  ice  layer  shape  and 
appearance  remained  similar  to  that  observed  under  steady-state  conditions; 
i.e.,  the  surface  was  smooth  and  level  across  the  channel.  The  end  of 
the  transient  was  indicated  by  reestablishment  of  steady  or  quasi-steady 
conditions  with  respect  to  ice  layer  thickness  and  cooling  plate 


! 
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temperature. 


The  variations  of  T and  T,  witli  time  were  monitored  by  the 

w b 

strip  chart  recorder.  T , the  bulk  water  temperature,  varied  only 

slightly  during  a given  run  while  the  wall  temperature,  T , usually 

w 

varied  30°C  to  40°C.  The  ice  layer  thickness,  S,  was  obtained  by  a 
mechanical  probe  technique  somewhat  modified  from  that  used  in  the 
steady-state  experiments.  The  probe  tip  was  positioned  approximately 
0.3  mm  to  0.6  mm  below  the  steady-state  ice  layer  surface  at  the  start  of 
a run.  When  the  cooling  plate  temperature  transient  began,  the  ice  layer 
grew  and  the  layer  surface  approached  the  tip  of  the  probe.  When  it  was 
determined  (by  observation  through  the  magnifying  lens)  that  the  layer 
had  contacted  the  tip,  the  probe  was  quickly  moved  downward  0.127  mm 
(0.005  inch).  This  process  was  repeated  as  the  layer  grew  and  approached 
the  probe  at  each  new  position.  Simultaneously,  the  thermocouple  output 
of  the  probe  was  recorded  on  a strip  chart  recorder.  In  this  way,  the 
chart  trace  indicated  a sharp  probe  tip  temperature  rise  each  time  the 
probe  was  pulled  downward  out  of  the  thermal  boundary  layer.  The  trace 
of  probe  tip  temperature  represented  a record  of  S,  ice  thickness,  vs. 
time,  in  that  the  time  elapsed  between  successive  discontinuities  in  the 
strip  chart  trace  indicated  the  time  required  for  the  layer  to  grow 
0.127  nmi.  Normally,  this  procedure  yielded  from  10  to  20  thickness 
measurements  in  a single  transient  run,  depending  upon  how  much  growth 
the  layer  experienced. 

2 . Combination  Tlnto  Temperature  an d Interface  Heat  Flux  Transients 
Runs  of  this  type  differed  from  the  plate  temperature  transient 
runs  in  that  the  layer  growth  was  induced  by  changes  in  both  interface 
convective  heat  flux  and  cooling  plate  temperature.  The  procedure  used 
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was  to  build  a steady-state  ice  layer  on  the  cooling  plate  with  the 
coolant  bin  fully  loaded  with  dry  ice.  The  interface  heat  flux  was 
maintained  at  a high  level  by  maintaining  high  volumetric  flow  rates 
through  the  channel,  so  that  the  steady-state  ice  layer  thickness  would 
be  small.  In  the  combination  plate  temperature  and  interface  heat  flux 
transient  runs,  the  initial  ice  thickness,  S^,  was  measured  directly  by 
the  mechanical  probe  method.  After  this  measurement  had  been  made  (time 
t = 0) , the  water  flow  rate  was  decreased  to  some  intermediate  value  by 
rapidly  closing  the  water  throttling  valve.  This  in  turn  decreased  the 


interface  heat  flux  and  initiated  ice  layer  growth.  The  change  in  heat 
flux  occurred  over  such  a short  period  of  time  (about  1 second)  that, 
for  all  practical  purposes,  it  was  a step  change.  For  times  t > 0,  the 
flow  rate  remained  constant,  as  did  the  interface  heat  flux,  except  for 
a slight  increase  due  to  the  small  variation  in  bulk  water  temperature 
which  occurred  during  the  runs. 

The  step  decrease  in  interface  heat  flux  at  time  t = 0 always 
caused  the  cooling  plate  surface  temperature  to  drop.  High  heat  fluxes 
tended  to  drive  the  plate  temperature  higher  while  the  plate  temperature 
dropped  under  conditions  of  low  interface  heat  flux.  This  behavior  could 
be  caused  by  a decrease  in  film  coefficient  in  the  coolant  bin  due  to 
increased  dry  ice  sublimation  (and  thus  increased  voidage)  at  the  upper 
surface  of  the  cooling  plate  under  high  interface  heat  flux  conditions. 
Because  of  this  dependence  of  plate  surface  temperature  on  interface 
heat  flux,  ice  layer  growth  transients  initiated  only  by  a change  in  heat 
flux  were  not  possible.  Therefore,  ice  layer  growth  was  due  to  a step 
decrease  in  Interface  heat  flux  at  t = 0 followed  by  a slower  drop  in 
cooling  plate  surface  temperature.  Plate  temperature  and  bulk  water 


A 


73 


temperature  throughout  these  runs  were  recorded  on  strip  chart  recorders. 
Ice  layer  thickness  measurements  were  obtained  by  the  same  method  as 
that  described  for  the  plate  transient  runs.  Table  5 (Chapter  IV) 
summarizes  the  transient  experiments. 


3.  Uncertainty  in  the  Transient  Measurements 

The  uncertainties  involved  in  the  measurements  of  water  flow  rate, 
bulk,  temperature,  and  plate  temperature  were  stated  in  the  description  of 
these  measurements.  The  estimated  uncertainties  also  apply  to  the  raw 
flow  rate  and  temperature  measurements  obtained  in  the  transient  runs. 
However,  the  uncertainty  in  these  raw  measurements  does  not  completely 
reflect  the  total  uncertainty  involved  in  the  inputs  to  the  finite  differ- 
ence model  such  as  St  and  q*  . The  additional  complications  associated 
with  the  estimation  of  these  parameters  will  be  discussed  later.  The 
uncertainty  in  the  raw  measurement  of  S,  the  ice  layer  thickness,  during 
the  transient  runs  deviated  from  the  estimated  value  of  +0.05  mm  stated 
in  the  section  on  ice  thickness  measurement.  The  +.05  mm  applies  primarily 
to  the.  steady-state  measurement  of  S . In  the  technique  for  measuring 
transient  ice  layer  thickness,  it  was  necessary  that  the  probe  remain 
near  to  the  ice  layer  interface  for  an  extended  period  of  time  (the  probe 
tip  was  never  more  than  0.127  mm  from  the  ice  surface  during  most  of  the 
transient)  and  some  local  melting  of  the  ice  layer  was  observed  in  the 
immediate  vicinity  of  the  probe.  Figure  18,  a and  b,  shows  schematically 
the  type  of  layer  deformation  which  occurred.  In  this  figure,  S'  denotes 
the  distorted  layer  thickness,  while  S is  the  thickness  of  the  undisturbed 
layer  under  the  same  conditions.  Thus,  it  is  apparent  that  the  thickness 
measurements  made  in  the  transient  experiments  indicate  S'  rather  than 
the  true,  undisturbed  thickness,  S . A simple  technique  was  therefore 
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Figure  18.  Sketch  of  Ice  Layer  Distortion  During 
Transient  Thickness  Measurement 
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developed  which  would  predict  S,  the  true  ice  thickness,  given  the 
measured  thickness  S'  and  the  flow  conditions  in  the  channel. 

The  phenomenon  described  above  is  thought  to  arise  due  to  a dis- 
torted water  flow  pattern  near  the  probe  tip  which  in  turn  causes  a 
local  increase  in  the  convective  film  coefficient  at  the  interface.  If 
this  hypothesis  is  adopted,  the  interface  film  coefficient  during 
transient  ice  thickness  measurement  becomes: 


h’  = h + Ah  , 


in  which  h is  the  film  coefficient  for  the  given  flow  conditions  over 
the  undisturbed  layer  and  Ah  is  the  incremental  increase  in  film 
coefficient  due  to  the  probe  induced  flow  disturbance.  During  an  ice 
layer  growth  transient,  the  conditions  in  the  layer  and  at  the  ice  layer 
surface  are  approximately  as  shown  in  Figure  19.  An  energy  balance  at 
the  layer  interface  results  in  the  following  relationships. 

For  the  undisturbed  portion  of  the  interface. 


„ II  _ -n  uj 

" qc  ~ dF  * 


and  for  the  disturbed  interface  immediately  above  the  probe  tip. 


. 0 dS ' 

qc  = pi  , 


in  which  q'  is  the  heat  flux  to  the  distorted  portion  of  the  layer 
c 

surface;  i.e.. 
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Figure  19 


Distorted  and  Undistorted  Interface 
During  Ice  Layer  Growth  Transient 
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q’  = h'  (T  -T  ) = (h+Ah) (T  -T  ) 
c b m b m 


(41) 


In  order  to  develop  an  approximate  expression  for  S in  terms  of  Ah 
and  S',  the  following  assumptions  were  made: 


1.  The  temperature  profile  in  the  layer  is  approximately  linear. 


thus ; 


9T 

9y 


(T  -T  ) 
m w 


(42) 


y=S 


and 


8T 

3y 


(T  -T  ) 
m w 


S' 


(43) 


y=S' 


This  assumption  is  certainly  satisfactory  for  the  order  of  approximation 


sought  here  because  of  the  small  ratio  of  specific  to  latent  heat 

„-l 


c/fL  = 0. 006°C 


dS  dS  ^ 

2.  The  interface  velocities  ^ and  are  nearly  equal: 


dS  = dSl 
dt  dt 


(44) 


which  follows  from  the  following  reasoning.  If  the  two  interface 
velocities  differed  significantly  in  magnitude,  the  distortion  of  the 
interface  would  become  more  pronounced  as  the  layer  growth  proceeded; 
this  was  not  observed  in  the  experiments. 
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Using  these  assumptions  and  combining  Equations  (34)  and  (40)  yields: 


(T  -T  ) 


(T  -T  ) 


_ n"  = V;  ^ 

qc  S' 


(45) 


By  noting  that. 


(h+Ah) (T  -T  ) 


b m 


(46) 


and 


h(T  -T  ) * 

b m 


(47) 


Equation  (45)  becomes: 


Si  = rS'  Ah,  (Tb  Tm) 
S 1 - ( k ; (T  -T  ) 


(48) 


Equation  (48)  relates  S,  the  true,  undisturbed  ice  layer  thickness  to  the 
measured  thickness  S'  and  the  incremental  increase  in  film  coefficient 
Ah  . At  a given  time  during  a transient  run,  all  quantities  on  the  right- 
hand  side  of  Equation  (48),  except  Ah,  are  known.  Ah  was  determined 
through  a series  of  steady-state  experiments.  Under  steady-state  conditions, 


h(T  -T  ) = k 
b m 


(T  -T  ) 
m w 


(for  the  undisturbed  interface)  (49) 


and 


(T  -T  ) 

(h+Ah) (T  -T  ) = k ra  w 


b m 


s' 


(for  the  disturbed  interface) . (50) 


■ — 


h 
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Thus,  for  the  steady-state  ice  layer. 


h+Ah  S 
h = s'  • 


(51) 


Both  S and  S'  were  determined  by  building  a steady-state  ice  layer  on 
the  cooling  plate.  The  probe  was  quickly  advanced  toward  the  layer  surface 
and  the  thickness  S of  the  undisturbed  layer  noted.  The  probe  was  allowed 
to  remain  in  this  position  for  10  to  15  seconds  during  which  time  the  layer 
receded  slightly  due  to  Ah,  the  local  increase  in  film  coefficient  caused 
by  the  probe.  Then,  the  probe  was  again  moved  upward  until  contacting  the 
interface  and  S'  was  measured.  Ah  was  then  computed  from  Equation  (51). 
Ah  was  found  to  be  primarily  a function  of  volumetric  flow  rate  in  the 
channel;  therefore,  measurements  of  Ah  were  made  at  the  three  volumetric 
flow  rates  used  in  the  transient  runs.  The  quantity  (h+Ah)/h  is  given 
in  Table  2 for  these  three  flow  rates. 


TABLE  2 


Tabulation  of  the  Quantity  (h+Ah)/h 


sec 

h+Ah 

h 

-3 

0.497  X 10 

1.008 

-3 

0.639  X 10 

1.022 

0.860  X 10-3 

1.038 

j 2,  j A * ' »■. 


CHAPTER  IV 


RESULTS  AND  DISCUSSION 

A.  Reduction  of  Experimental  Data 

The  results  of  the  steady-state  heat  flux  measurements  are  used  to 
obtain  a correlation  for  interface  film  coefficient  as  a function  of 
flow  conditions  within  the  test  channel.  Also,  the  analytic  representa- 
tic  of  the  wall  temperature  variations  obtained  in  the  transient  runs 
are  discussed  in  some  detail. 

1.  The  Correlation  for  Interface  Film  Coefficient 
The  iriterfacial  film  coefficient,  h,  was  assumed  to  be  dependent 
upon  only  two  flow  variables,  the  volumetric  rate  of  water  flow  through 
the  test  section  and  the  bulk  water  temperature  entering  the  test  section. 
Curves  were  generated  for  the  film  coefficient  as  a function  of  bulk 
water  temperature  at  constant  flow  rate  for  three  different  values  of  the 
flow  rate  Q,  and  bulk  temperatures  between  12°C  and  24°C.  A dimensional 
correlation  of  the  form, 

h = h[T,  ] , (constant  Q)  (52) 

b 

was  found  to  be  the  most  convenient  for  use  in  the  transient  experiments 
because  the  volumetric  flow  was  always  held  constant  in  a given  transient 
run.  Dimensionless  correlations  of  the  form, 


6 r\> 
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Nu  = Nu[Re,  Pr]  (53) 

were  found  to  be  unsatisfactory  because  at  constant  flow  rate,  both  Re 
and  Pr  vary  with  bulk  temperature  due  to  the  sizeable  viscosity 
variation  of  water  in  the  range  12°C  to  24°C.  The  three  volumetric  flow 
rates  and  the  range  of  bulk  water  temperatures  used  in  the  correlation 
coincide  with  the  flow  rates  and  temperatures  used  in  the  transient  runs. 

The  results  of  the  steady-state  heat  flux  experiments  are  shown  in 
Figures  20,  21,  and  22  for  each  of  the  three  water  flow  rates.  In  these 
figures,  the  uncertainty  associated  with  each  of  the  discrete  data  points 
depends  upon  the  cooling  plate  temperature,  the  bulk  water  temperature, 
and  the  steady-state  ice  layer  thickness  at  the  time  of  each  heat  flux 
measurement  and,  therefore,  varies  from  data  point  to  data  point.  The 
maximum  value  of  the  uncertainty  in  h occurring  in  the  measurements  at 
each  flow  rate  is  given  in  the  figures. 

The  solid  curves  in  Figures  20,  21,  and  22  represent  the  results 
of  the  dimensionless  correlation, 

Nu  = 0.0155  /Pr  Re0,83  , (54) 

an  equation  given  by  Kays  (19)  for  fully  developed,  turbulent  flow  in 
cylindrical  tubes.  In  Equation  (54),  the  hydraulic  diameter  of  the  flow 
channel  was  used  in  forming  the  dimensionless  parameters  Nu  and  Re  . 

This  correlation  gives  a rough  estimate  of  the  film  coefficient  at  the  ice- 
water  interface,  useful  for  comparison  with  the  results  of  the  experimental 


correlation. 


INTERFACE  FILM  COCFFICIEJ 


EQUATION  (54) 

K'P:  11,140  14,250 
Pr  6.4  - 8.0 
Q = 0.49/  x 10'3  m3/sec 
V = 0.348  ni/sec 


MAXIMUM  UNCERTAINTY:  + 4.8% 

I L 1 I . „..l I I 1 I _ I 

4 16  18  20  22  24 

BULK  TEMPERATURE  (deg  Cl 


Figure  20. 


Interface'  Film  Coefficient  Versus  Bulk 
Water  Temperature,  Q * 0.407  X 10''  mV 


EQUATION  (54) 

RC;  18,750  - 23.450 
Pr.  6.70  - 8.50 
Q = 0.860  x io  ^ m^/sec  " 
V = 0.602  m/sec 

MAXIMUM  UNCERTAINTY.  + <1.7% 

1 1 1 i 1 i i . i i 

M 16  IS  20  22 

BULK  TEMPERATURE  (dixj  C) 


Interface  Film  Coofflcienl  Versus  Bulk 
Water  Temperature,  Q » 0.860  X 10“-*  ni-^/ 
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The  results  shown  in  Figures  20  through  22  indicate  that  for  each 
flow  rate,  the  flow  at  the  ice  thickness  measurement  position  is  not 
fully  developed  thermally.  This  is  not  unexpected  because  cooling  takes 
place  for  a length  of  only  3.5  hydraulic  diameters  before  the  point  of 
heat  flux  measurement.  Also,  it  should  be  pointed  out  that  the  commonly 
used  dimensionless  fluid  property  and  flow  parameters  Pr  and  Re  behave 
approximately  as 

RePr  = CONSTANT  (55) 

in  each  of  Figures  20,  21,  and  22.  Therefore,  when  moving  to  the  left 

on  the  abscissa  (decreasing  T.  ),  Pr  is  increasing  while  Re  is 

b 

decreasing  according  to  Equation  (55).  The  experimental  value  of  h 
decreases  with  decreasing  bulk  temperature  for  basically  two  reasons. 

First,  Re  in  the  channel  is  decreasing  as  bulk  temperature  is  lowered. 
This  should  cause  the  film  coefficient  to  decrease  even  though  Pr  is 
increasing,  because  Nu  is  generally  a stronger  function  of  Re  than 
Pr  . Secondly,  because  Pr  is  increasing  with  decreased  bulk  temperature, 
the  temperature  profile  should  be  more  nearly  fully  developed  at  lower 
temperature,  resulting  in  decreased  film  coefficients.  The  effect  of 
increased  Pr  on  the  rate  of  thermal  development  in  turbulent  channel 
flow  is  demonstrated  by  Kays  (19).  The  downward  slope  of  the  data  for 
increased  Pr  (decreasing  T.  ) is  seen  to  be  less  and  less  apparent  as 

D 

flow  rate  is  increased;  in  fact,  in  Figure  22,  the  effect  is  nearly 
undetectable.  This  indicates  that  the  effects  of  Pr  on  the  flow 
development  may  become  less  pronounced  at  higher  Re  (higher  flow  rate). 


, jl 
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For  the  purpose  of  comparison,  the  data  of  Figures  20,  21,  and  22 
are  replotted  in  Figure  23.  In  this  figure,  the  solid  lines  represent 
least  square  (polynomial)  curve  fits  to  the  three  data  sets.  The 
absolute  spread  in  the  data  is  seen  to  be  largest  for  the  top  (high  flow 
rate)  curve  in  Figure  23.  This  reflects  the  difficulties  in  maintaining 
sufficiently  low  cooling  plate  temperatures  and  adequate  ice  layer  thick- 
ness in  the  presence  of  the  high  interfacial  heat  fluxes.  Table  3 
provides  a brief  summary  of  the  results  of  the  steady-state  heat  flux 
measurements.  In  this  table,  "V"  is  the  average  flow  velocity  in  the 
flow  channel  just  upstream  of  the  cooling  plate  leading  edge,  (DEV) rmg 
is  the  root  mean  square  deviation  of  the  data  from  the  curve  generated 

by  least  squares,  and  (DEV)  is  the  maximum  percentage  deviation  from 

mdx 

the  least  squares  curve. 

The  results  tabulated  in  Table  3 show  that  the  prediction  of 
interfacial  film  coefficient  is  probably  not  in  error  by  more  than  +5%. 

2 . Cooling  Plate  Temperature  Variations 

During  the  ten  transient  ice  growth  experiments,  cooling  plate 

surface  temperatures  were  measured  as  a function  of  time.  It  was  necessary 

that  these  temperature  variations  be  placed  in  a form,  T (t) , which 

w 

would  permit  reliable, accurate  wall  temperatures  (and  derivatives)  to  be 
computed  for  any  given  time,  t . 

Due  to  the  rapid  plate  surface  temperature  drop  during  the  first 
few  moments  of  the  transient  runs  and  the  asymptotic  approach  toward  a 
steady-state  plate  temperature,  polynomial  curve  fits  by  the  method  of 
least  squares  of  the  form 


BULK  TEMPERATURE.  °C 


Figure  23.  Results  of  the  Interface  Film 
Coe f f i c l en t Co r re  1 a t i on 
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did  not  satisfactorily  represent  the  experimental  data.  Therefore,  the 
plate  temperature  data  was  represented  by  rational  functions  of  the  form 


T 

w 


Aj+A,t+A  t2 
1+B1t+B^t2. 


.+A 


J+l 


,+B  t 

K 


j 

t 


k 

1+  E B,  tK 
k-1  k 


where  the 


and  B,  are  constants  and  J and  K 
k 


1 and  7. 


(57) 


each  ranged  between 


This  form  of  equation  allows  much  more  flexibility  than  Equation 
(56)  in  the  type  of  functional  relationship  which  is  to  be  represented. 

This  is  because  a simple  polynomial  expression  such  as  Equation  (56) 
cannot  have  any  poles  and  thus,  cannot  easily  approximate  functions  with 
singularities.  Equation  (57),  on  the  other  hand,  may  have  as  many  as  K 
distinct  poles.  This  characteristic  makes  the  rational  function  ideal 
for  the  representation  of  many  critical  phenomena  [see  Baker  (21)],  and 
functions  which  are  singular  or  nearly  singular.  However,  a problem 
arises  in  that  the  poles  of  Equation  (57)  may  occur  anywhere  within  the 
domain  of  the  original  function,  and  these  singularities  may  be 
undesirable. 

The  constants  A.  and  B,  were  chosen  so  that  Equation  (57) 

J k 

would  collocate  or  coincide  with  the  experimental  data  points  at  J + K + 1 
points,  spaced  evenly  with  respect  to  the  temperature  measurements.  For 
a given  transient  run,  these  coefficients  were  computed  for  all  permuta- 
tions and  combinations  of  J and  K from  J » 1,  K “ 1,  to  J-7, 

K ■ 7 resulting  in  49  different  rational  function  correlations  for  T 

w 


: 


* x.i 
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From  this  group,  one  was  chosen  which  resulted  in  the  best  fit  of  the 
data  and  which  had  no  poles  within  the  time  domain  of  interest.  Figure 
24  shows  the  fit  obtained  for  transient  run  5,  in  which  the  rational 
function  used  had  a numerator  of  degree  1 and  a denominator  of  degree  4. 
The  fit  is  seen  to  be  good  for  all  values  of  t including  small  t where 
plate  temperature  is  dropping  rapidly. 


dSt 


dT 


The  derivative  requires  the  evaluation  of  from  the 

dt*  dt* 

transient  experimental  data,  this  term  was  evaluated  through  differenta- 


tion  of  Equation  (57): 

«.  - i 1 w 


j-i 


K . ] 

Z kB  tk_i 

Z A ....  t^ 

k-i  k 

>0  J+1 

dt  K 

1 + I B t 
k-i  k 


K l 
1 + Z B,  t 

k-1  k 


(58) 


The  results  of  the  rational  function  approximations  for  all  ten 
transient  runs  are  shown  in  Table  4.  In  this  table, 


RMS  ERROR  = 


I 1/2 


I [Twi-R(t.)]‘ 
i=l  1 


(59) 


in  which  T . is  the  measured  wall  temperature  at  t = t.,  R(t.)  is  the 

wi  ii 

rational  function  approximation 
points  used  in  the  correlation. 


rational  function  approximation  to  T . , and  I is  the  number  of  data 

wi 


B.  Analytical  Results 

1.  Results  of  the  Finite  Difference  Model 

Sample  results  of  the  finite  difference  calculation  for  the  solid 
layer  thickness  were  generated  for  four  cases.  It  was  assumed  in  all 
cases  that  the  solid  layer  is  in  a state  of  equilibrium  prior  to  the 
time  t = 0 . The  initial  temperature  profile  in  the  layer  then,  is 


Figure  24.  The  Rational  Function  Curve  Fit  for  T , Transient  Run 
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given  by: 


0Q(y*)  = y* 


(t*  = 0) 


(60) 


i.e.,  the  initial  temperature  profile  is  linear.  This  assumption, 
however,  need  not  be  made  as  the  model  may  allow  initial  conditions  other 
than  steady-state.  At  time  t*  = 0,  either  the  wall  temperature  or  the 
interface  heat  flux  changes  in  a stepwise  manner,  thus,  initiating 
solid  layer  solidification  or  melting.  The  initial  values  of  St  and 
q*  (before  layer  growth  or  decay  begins  , i.e.,  for  t*  < 0)  are 

indicated  by  the  subscript  "i",  while  the  unsubscripted  symbols  refer 
to  the  values  of  the  parameters  for  times  t*  _>  0 . Note  that  for  a 
steady-state  layer  initially, 


Sti  ‘ ’J 


(61) 


and  that  the  final  steady-state  value  of  S*  is  given  by, 


St 


S*  = — (at  steady-state) 


(62) 


In  Figures  25  and  26,  the  behavior  of  the  solid  layer  is  shown  in 
terms  of  the  dimensionless  coordinates  S*  and  t*  . In  these  figures 
the  interface  heat  flux  is  held  constant  while  the  wall  temperature  is 
decreased  promoting  solid  layer  growth;  while  in  Figure  26,  the  wall 
temperature  is  increased  causing  the  layer  to  melt.  Both  figures  clearly 
show  the  significant  effect  which  the  parameter  St  has  on  the  layer 
transient  behavior.  In  Figure  25  (solidification),  the  interface  velocity 
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dS*/dt*  is  initially  zero;  the  initial  velocity  appears  to  be  non-zero 
only  because  the  time  scale  in  Figure  25  is  too  coarse  for  the  initial 
period  of  rapid  acceleration  to  be  clearly  seen.  The  most  important 
aspect  of  the  curves  in  Figure  25  is  the  effect  of  St  on  the  time 
required  for  the  layer  to  assume  a new  steady-state  thickness.  This 
time  is  substantial ly  increased  for  increased  Stefan  number;  in  fact, 
even  the  basic  shape  of  the  curves  change  with  increased  St,  eventually 
assuming  the  characteristic  "square  root  shape"  of  the  Stefan  solution, 

(2),  for  very  large  values  of  St  . 

In  order  to  provide  some  insight  into  this  behavior,  the  dimension- 
less solid  layer  temperature  profiles  are  shown  in  Figures  27  and  28  for 
Stefan  numbers  of  2 and  10,  respectively,  and  conditions  corresponding 
to  those  in  Figure  25.  In  these  two  figures,  the  approach  toward  steady- 
state  is  indicated  by  the  decreasing  degree  of  curvature  in  the  temperature 
profiles  (for  increased  t*)  or,  more  specifically,  the  approach  of  the 
slope  at  y*/S*  **  1 (phase  boundary)  toward  unity.  The  reason  for  the 
rapid  attainment  of  steady-state  conditions  for  small  St  is  now  clear, 
for  the  temperature  profile  in  the  layer  and  its  approach  toward  the 
steady-state  is  significantly  affected  by  the  parameter  St  . Note  that 
in  Figure  27,  for  St  = 2,  the  temperature  profile  is  almost  linear  at 
t*  » 1,  while  in  Figure  28,  for  St  = 10,  the  temperature  profile  possesses 
considerable  curvature  even  for  t*  = 10  . 

In  Figure  26,  similar  behavior  is  exhibited  by  melting  solid  layers. 

The  initial  interface  velocity  is  again  zero  and  the  melting  rate  is  seen 
to  increase  with  decreased  St  . For  very  small  St,  the  dimensionless 
layer  thickness,  S*,  becomes  approximately  a linear  function  of  t*  as 
is  apparent  from  the  interface  energy  balance  condition  (5f) ; i.e.  , 

. 

“ — — ■ ■ — - 
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dS* 
dt  * 


-q* 


S*  a l-q*t* 


> St  > 0 


((.3) 


(64) 


It  should  also  be  noted  that  Equation  (64)  Is  the  limit  of  Equation  (10)  , 
the  asymptotic  solution,  as  St  0 . As  in  Figure  25,  the  time  required 
for  the  layer  thickness  to  reach  a new  steady-state  value  is  again 
largely  dependent  upon  St  . Steady-state  is  seen  to  be  attained  more 
quickly  for  small  St  than  for  large  St,  even  though  more  of  the  layer 
must  melt  when  St  is  small. 

In  Figures  29  and  30,  the  solid  layer  behavior  is  shown  for  the 
constant  wall  temperature  case  where  interface  heat  flux  is  decreased 
(Figure  29)  and  increased  (Figure  30)  resulting  in  solidifying  and  molting 
solid  layers,  respectively.  In  Figure  29,  the  interface  velocity  is 
non-zero  at  t*  = 0 . The  layer  growth  profiles  in  Figure  29  are  similar 
in  shape  to  those  in  Figure  25  for  constant  heat  flux;  however,  in  this 
case,  it  is  the  interfacial  heat  flux  parameter,  q*,  which  governs  the 
rate  at  which  steady-state  conditions  are  attained.  In  Figure  30,  the 
heat  flux  is  increased  at  t*  “ 0 so  that  the  solid  layer  melts  down  to 
some  fraction  of  its  original  thickness.  Once  again,  very  large  values 
of  the  heat  flux  cause  the  layer  to  melt  at  a rate  which  is  almost 
constant,  and  reach  a new  steady-state  very  quickly.  The  one  obvious 
conclusion  which  may  be  drawn  from  Figures  25,  26,  29  and  30  is  that  both 
St  and  q*  not  only  govern  the  magnitude  of  the  final  steady- state 
layer  thickness,  but  more  importantly,  these  parameters  govern  the 
character  of  the  solution  S*(t*)  and  the  duration  of  the  transient. 


In  order  to  illustrate  the  accuracy  and  effectiveness  of  the 
finite  difference  model,  at  least  for  a simple  case  of  constant  wall 


temperature  solidification,  the  numerical  solution  is  compared  with  the 
exact  solution  of  Stefan  (2).  An  inconsistency  arises  in  this  comparison 
due  to  the  fact  that  the  numerical  solution  as  formulated  in  Chapter  II 
assumes  that  the  solid  layer  is  always  of  finite  thickness,  while  in 
the  Stefan  solution,  the  solid  layer  is  of  zero  thickness  as  t = 0 . 

The  comparison  of  the  exact  and  numerical  solutions  is  therefore  made  by 
initially  computing  the  thickness  variation  from  the  exact  solution  and 
continuing  in  this  way  until  the  solid  layer  is  Sq  in  dimension,  at 
which  time  the  finite  difference  solution  is  begun.  In  terms  of  the 
pertinent  variables  of  the  present  analysis,  the  Stefan  solution  becomes: 


l_2/t*+l/4b 
erf (b) 


S*  = 4b  t*  + 1 , 


where  the  constant  b is  determined  from: 


,2  b^  ...  . St 
b e erf(b)  = — . 


For  convenience,  St  was  set  equal  to  0.5923  so  that  b = 1/2  ; then. 
Equation  (66)  becomes: 


S*  = /t*+l  . 
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The  results  of  Equation  (68)  and  the  finite  difference  solution 
are  compared  in  Figure  31.  This  figure  indicates  that,  at  least  in  this 
simple  case,  there  is  essentially  no  difference  between  the  exact  and 
numerical  solutions.  The  results  of  this  comparison  and  the  results  of 
others  to  follow  indicate  that  the  finite  difference  solution  technique 
does  indeed  converge  to  the  exact  solution  of  the  system.  Equations (5) . 

2 . Evaluation  of  the  Approximate  Solution  Techniques 

It  was  mentioned  previously  that  the  extent  to  which  the  approximate 
techniques  of  solution  (asymptotic  and  integral)  are  capable  of  represent- 
ing the  exact  solution  can  only  be  established  through  quantitative 
comparison  with  methods  which  are  believed  to  be  exact.  Below,  the 
overall  effectiveness  of  each  of  the  approximate  techniques  is  investigated. 

In  Figure  32,  the  asymptotic  solution  for  small  St  is  compared 
with  the  finite  difference  numerical  solution  for  Stefan  numbers  of  1/50, 
1/10,  and  1.  From  this  figure,  the  effect  of  St  on  the  accuracy  of  the 
asymptotic  solution  is  evident.  For  St  as  large  as  1/10,  the  approximate 
solution  represents  the  true  solution  with  little  error,  while  for  St  = 1 , 
the  accuracy  of  the  asymptotic  solution  rapidly  decays.  The  capability 
of  the  asymptotic  solution  in  yielding  accurate  results  depends  primarily 
upon  how  well  the  assumed,  linear  solid  layer  temperature  profile 
approximates  the  actual  temperature  profile.  The  actual  temperature 
profiles  (computed  from  the  full  finite  difference  numerical  solution) 
are  shown  in  Figure  33  for  the  conditions  of  Figure  32  and  Stefan  numbers 
of  l and  1/50.  The  curves  in  this  figure  clearly  indicate  that  for  any 
time,  t*,  the  profile  for  small  St  is  a substantially  better 
approximation  to  the  steady-state,  linear  profile  than  the  profile  for 
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Figure  33.  Dimensionless  Solid  Layer  Temperature 
Profiles,  St  = 1/30  and  1 
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large  St  and  it  therefore  follows  that  the  asymptotic  solution  yields 
more  accurate  results  when  St  is  small.  The  conclusion  drawn  from 
the  above  discussion  is  that  the  asymptotic  solution,  Equation  (10), 
does  yield  an  accurate,  simple  means  of  estimating  the  solid  layer 
thickness  variation  for  small  values  of  the  parameter  St  . Of  course, 
for  instances  in  which  St  and/or  q*  are  not  constant,  Equation  (9) 
must  be  used  in  computing  S*  rather  than  Equation  (10).  The  practical 
utility  of  Equations  (9)  and  (10)  depend  primarily  upon  the  specific 
substance  and  wall  temperature  variation  in  the  problem  under  considera- 
tion i.e.  the  magnitude  of  St  = c(T  -T  )/2  . Ice  (&/c  2 160°C)  is 

m w 

one  substance  for  which  the  asymptotic  solution  could  prove  very  useful 

provided  that  the  subcooling  (T  -T  ) in  the  ice  layer  is  small 

in  w 

compared  with  H/c  ; this  is  often  a good  assumption. 

In  Figures  34  through  36,  the  approximate  integral  solution, 
Equation  (19) , is  shown  plotted  along  with  the  exact  numerical  solution 
for  Stefan  numbers  2,  5,  and  10.  Unlike  the  asymptotic  solution,  the 
integral  solution  tends  to  yield  reasonably  accurate  solid  layer 
thickness  results  independent  of  St  . In  all  of  the  plots  the  maximum 
error  in  the  integral  solution  tends  to  be  about  5%,  the  larger  errors 
occurring  for  small  time  t*  . Evidently,  the  parabolic  temperature 
profile  assumed  in  Equation  (18)  does  reasonably  well  in  approximating 
the  actual  solid  layer  temperature  profile  and  does  especially  well  for 
large  time,  t*,  when  the  profile  tends  toward  the  linear,  steady-state 
shape.  From  these  results  it  appears  that  the  integral  solution  tech- 
nique may  be  as  general  a solution  method  and  nearly  as  accurate  as  the 
finite  difference  solution  and  this  is  indeed  true  for  many  situations. 
However,  the  form  of  Equation  (19)  is  such  that  dSt/dt*  and  dq*/dt* 


are  required  for  Che  integral  solution  of  problems  in  which  St  and  q* 
are  both  time  dependent.  This  factor  can  make  the  integral  technique 
somewhat  inconvenient  when  neither  wall  temperature  nor  interface  heat 
flux  are  known  precisely  as  functions  of  time;  in  such  situations  the 
finite  difference  solution  is  the  more  appropriate. 

C.  Evaluation  of  the  Mathematical  Model 

In  the  previous  sections  the  basic  methods  of  solution  (asymptotic, 
integral,  and  finite  difference)  were  tested  and  evaluated  with  respect 
to  their  ability  to  generate  accurate  solutions  to  the  general  mathematical 
model  as  formulated  in  Chapter  II.  In  this  section,  the  model  itself  and 
the  assumptions  and  simplifications  made  in  its  development  are  tested 
and  evaluated  via  the  results  of  the  transient  ice  growth  measurements 
described  in  Chapter  III. 

The  results  of  the  ten  transient  tests  are  compared  with  the  pre- 
dictions of  the  finite  difference  computer  model  in  Figures  37  through  46, 
while  Table  5 indicates  the  type  of  transient  and  the  conditions  present 
in  the  flow  channel  for  each  of  the  ten  runs.  The  computer  predictions 
are  indicated  by  the  solid  curves  while  the  results  of  the  transient  ice 
thickness  measurements  are  represented  by  the  data  points.  The  range  of 
uncertainty  shown  with  each  data  point  corresponds  to  the  +0.05  mm 
uncertainty  in  the  transient  ice  thickness  measurement  after  the 
correction  was  made  for  ice  layer  distortion.  In  Figures  31  through  41 
(plate  transient  runs) , the  agreement  between  the  experimental  and 
computer  results  (the  computer  predictions  are  indicated  by  the  solid 
curves)  is  seen  to  be  good  for  all  values  of  the  dimensionless  time  t*  . 
There  is  some  tendency  for  the  model  to  overpredict  the  ice  layer  growth 


TABLE  5 


Summary  of  Conditions  for  Transient  Runs 


Type  of 
Transient 

Flow  Rate 

x3iooo 

(m  /sec) 

Bulk  Temp. 
(°C) 

PLATE 

0.639 

17.6 

PLATE 

0.639 

19.0 

PLATE 

0.639 

18.3 

PLATE 

0.860 

15.7 

PLATE 

0.497 

20.  7 

COMB. 

0.639 

19.2 

COMB. 

0.639 

15.5 

COMB. 

0.497 

18.2 

COMB. 

0.639 

17.8 

COMB. 

0.860 

15.8 

Re  = VD.  /V 


'PLATE  denotes  plate  temperature  transient. 

'COMB."  denotes  combination  plate  temperature  and  inter 
face  heat  flux  transient. 


I 


Figure  44.  Results  of  Transient  Run  8 


rate  at  small  t*  while  underpredicting  it  at  large  t*  . The  same  is 
nominally  true  of  the  results  for  the  combination  plate  temperature  and 
interface  heat  flux  transient  tests  (Figures  42  through  46).  This  effect 
may  be  due  to  the  fact  that  constant  transport  properties  were  used  for 
the  ice  layer  in  the  model.  In  reality,  the  mean  thermal  conductivity 
in  the  layer  increased  by  as  much  as  15%  between  the  start  and  end  of  a 
transient  run,  while  the  conductivity  of  the  ice  near  the  plate  surface 
increased  by  as  much  as  30%  in  a run.  This  variation  would  cause  the 
actual  layer  growth  to  be  slower  than  predicted  at  the  beginning  of  the 
transient  and  more  rapid  toward  the  end.  The  value  for  ice  thermal 
conductivity  used  throughout  all  transient  runs  was  k = 2.44W/m-°C 
which  corresponds  to  a mean  ice  layer  temperature  of  -17.5°C  (20). 

The  conductivity  (and  hence,  diffusivity)  variation  may  also  be  partially 
responsible  for  the  fact  that  generally  better  agreement  was  obtained  in 
the  plate  transient  runs  because  the  time  averaged  mean  temperature  in 
the  layer  was  substantially  lower  than  -17.5°C  in  the  combination 
transient  runs.  It  should  also  be  pointed  out  that  the  interface  heat 
flux  correlation  used  in  the  formation  of  q*  introduces  an  uncertainty 
of  about  +5%  into  the  estimation  of  this  parameter.  The  only  other 
factor  which  may  have  made  some  contribution  to  the  deviation  between 
experiment  and  theory  was  the  analytic  representation  (and  subsequent 
differentiation)  of  the  wall  or  plate  temperature  variation  during  the 
transient  runs.  However,  due  to  the  overall  good  fit  yielded  by  the 
rational  function  approximations  and  to  the  smooth  derivatives  resulting 
from  their  differentiation,  it  is  believed  that  the  contribution  to  the 
deviations  seen  in  Figures  37  through  46  from  this  source  are  small. 
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D.  Final  Conclusions 

The  results  of  the  transient  experiments  and  the  comparison  with 
the  model  predictions  essentially  indicates  only  one  major  deficiency 
of  the  mathematical  model;  i.e.,  the  inability  to  take  into  consideration 
temperature  dependent  solid  layer  properties.  The  other  factors  mentioned 
above  which  may  have  contributed  to  the  discrepancy  between  experiment 
and  theory  result  from  uncertainties  in  the  model  input  data  rather  than 
from  any  of  the  fundamental  assumptions  made  in  the  model  formulation. 

In  an  attempt  to  estimate  the  seriousness  of  this  deficiency  with  respect 
to  the  practical  application  of  the  model,  one  factor  must  be  kept  in 
mind;  the  substance  used  in  the  present  experimental  investigation  is  not 
necessarily  a typical  one  with  respect  to  the  degree  to  which  its 
transport  properties  depend  on  temperature  near  the  fusion  temperature. 

It  is  true  that  dielectric  solids  such  as  ice  may  exhibit  a substantial 
variation  in  thermal  conductivity  with  temperature  [see  Dillard  (22)]. 
However,  most  substances,  pure  metals,  for  example,  which  commonly  occur 
in  melting  and  solidification  phenomena  have  thermal  conductivities  and 
dif fusivities  which  are.  not  nearly  as  dependent  on  temperature  as  are  the 
corresponding  properties  of  ice  [see  (23)  and  (24)].  It  is  therefore 
expected  that  the  failure  of  the  model  to  account  for  varying  transport 
properties  would  be  found  much  less  significant  in  experiments  conducted 
with  these  substances.  An  extension  of  the  present  work  might  attempt  to 
account  for  the  variation  of  the  solid  phase  properties  with  temperature 
in  the  finite  difference  model. 


In  the  final  evaluation  of  the  mathematical  model,  it  is  noted  that 
even  if  the  property  variation  of  ice  is  considered  typical,  the  results 
in  Figures  37  through  46  indicate  basically  good  agreement  between 


experiment  and  theory,  and  the  remaining  assumptions  and  simplifications 
made  in  the  model  formulation  (other  than  that  of  constant  solid  layer 
properties)  seem  to  be  justified  by  the  experimental  results.  It  is 
therefore  concluded  that  the  mathematical  model  and  calculation  procedure 
outlined  here  does  provide  an  effective  means  of  predicting  the  behavior 
of  the  simple,  two-phase  system  shown  in  Figure  1 under  conditions  of 
practical  engineering  interest. 


126 


REFERENCES 


1.  Cheung,  F.  B.  and  L.  Baker  Jr. , "Transient  Freezing  of  Liquids  in 
Tube  Flow,"  Nuclear  Science  and  Engineering.  Vol.  60,  p.  1-9,  1976. 

2.  Stefan,  J.,  "Uber  die  Theorie  der  Eisbildung,  insbesondere  iiber  die 
Eisbildung  im  Polarmeere,"  Annalen  der  Physik  und  Chemie,  Vol.  62, 
p.  269,  1891. 

3.  Portnov,  I.  C. , "Exact  Solution  of  the  Freezing  Problem  With 
Arbitrary  Temperature  Variation  on  the  Fixed  Boundary,"  Soviet 
Physics,'  Vol.  7,  No.  3,  p.  186,  1962. 

4.  Stephan,  K. , "The  Influence  of  Heat  Transfer  on  Melting  and  Solidifi- 
cation in  Forced  Flow,"  Int.  J.  Heat  Mass  Transfer,  Vol.  12,  p.  199-214, 
1969. 

5.  Goodman,  T.  R. , "The  Heat  Balance  Integral  and  Its  Application  to 
Problems  Involving  a Change  of  Phase,"  Trans.  A.S.M.E.,  Vol.  80, 
p.  335-342,  1958. 

6.  Libby,  P.  A.  and  S.  Chen,  "The  Growth  of  a Deposited  Layer  on  a Cold 
Surface,"  Int.  J.  Heat  Mass  Transfer,  Vol.  8,  p.  395-402,  1965. 

7.  Savino,  J.  M.  and  R.  Siegel,  "Experimental  and  Analytical  Study  of 
the  Transient  Solidification  of  a Warm  Liquid  Flowing  Over  a Chilled 
Flat  Plate ,"  NASA  TN  D-4015,  1967. 

8.  Biot,  M.  A.,  "New  Methods  in  Heat  Flow  Analysis  With  Application  to 
Flight  Structures,"  J.  Aeronaut.  Sci. , Vol.  24,  p.  857-873,  1957. 

9.  Lapadula,  C. , and  W.  K.  Mueller,  "Heat  Conduction  With  Solidification 
and  a Convective  Boundary  Condition  at  the  Freezing  Front,"  Int.  J. 

Heat  Mass  Transfer,  Vol.  9,  p.  702-704,  1966. 

10.  Bllenas,  J.  A.,  and  L.  Jiji,  "Variational  Solution  of  Axisynunetric 
Fluid  Flow  in  Tubes  With  Surface  Solidification,"  Journal  of  the 
Franklin  Institute,  Vol.  289,  No.  4,  April,  1970. 

11.  Springer,  G.  S. , and  D.  R.  Olson,  "Method  of  Solution  of  Axisymmet ric 
Solidification  and  Melting  Problems,"  Contributed  by  the  Heat  Transfer 
Division  for  presentation  at  the  Winter  Annual  Meeting,  Nov.  25-30, 

1962,  of  the  ASME. 

12.  Murray,  W.  D.  and  F.  Landis,  "Numerical  and  Machine  Solutions  of 
Transient  Heat  Conduction  Problems  Involving  Melting  or  Freezing," 

J.  of  Heat  Transfer,  ASME  Trans.  Series  C,  Vol.  81,  p.  106-112,  1959. 


127 


13.  E.  Friedman,  "An  Iterative  Procedure  for  Including  Phase  Change  in 
Transient  Conduction  Programs  and  its  Incorporation  into  the  Finite 
Element  Method,"  Submitted  through,  AIChE  HT  and  EC  Division  for  the 
1977  National  Heat  Transfer  Conference,  Solidification  and  Melting 
Heat  Transfer,  p.  182. 

14.  Schlicting,  Hermann,  Boundary  Layer  Theory,  McGraw-Hill,  New  York, 

1968. 

15.  Lapidus,  Leon  and  J.  H.  Seinfeld,  Numerical  Solution  of  Ordinary 
Differential  Equations,  Academic  Press,  New  York,  1971. 

16.  ftzisik,  M.  N. , Boundary  Value  Problems  of  Heat  Conduction,  International 
Textbook  Company,  Scranton,  PA.,  1968. 

17.  Beckwith,  T.  G. , and  N.  L.  Buck,  Mechanical  Measurements,  Addison- 
Wesley,  Menlo  Park,  California,  1973. 

18.  Krieth,  F. , Principles  of  Heat  Transfer,  Intext  Publishing  Company, 

New  York,  1976. 

19.  Kays,  W.  M. , Convective  Heat  and  Mass  Transfer,  McGraw-Hill,  New  York, 
1976. 

20.  Ratcliffe,  E.  H. , "The  Thermal  Conductivity  of  Ice-New  Data  on  the 
Temperature  Coefficient,"  Philosophical  Magazine,  Vol.  7,  No.  79, 
p.  1197-1203,  1962. 

21.  Saff,  E.  B.  and  R.  S.  Varga,  Pade  and  Rational  Approximation  (Theory 
and  Applications),  Academic  Press,  New  York,  1977. 

22.  Dillard,  D.  S. , and  K.  D.  Timmerhaus,  "Low  Temperature  Thermal 
Conductivity  of  Selected  Dielectric  Crystalline  Solids,"  Thermal 
Conductivity  Proceedings  of  the  8th  Conference,  Purdue  Univ. , 

West  Lafayette,  Indiana,  Oct.  7-10,  1968. 

23.  Ho,  C.  V.,  Powell,  R.  W. , and  K.  Y.  Wu,  "Thermal  Diffusivity  of  the 
Elements,"  Thermal  Conductivity  Proceedings  of  the  8th  Conference, 

Purdue  Univ.,  West  Lafayette,  Indiana,  Oct.  7-10,  1968. 

24.  Missenard,  Andre,  Conductivity  Thermique  Des  Solidcs,  Liquidcs,  Gaz  Et 
De  Leurs  Melanges,  Editions  Eyrolles,  Paris,  1965. 


i 


128 


APPENDIX  A 


FINITE  DIFFERENCE  APPROXIMATIONS 


Gr+1  _ flr 

l£*  ' t * ' ° + 0(At*> 


2 er  + 0r  - 20r 

-n_1-  n+1  — — + 0(Ay*2) 


3y*2 

Ay*2 

30  , 

0r+1  - 0r  , 

_ n+1  n-1 

3y* 

2Ay* 

30 

. ^-2 

3y* 

y*=S* 

dS* 

s*1^1  - s*r 

dt* 

At* 

+ 0(Ay*  ) 


4»Nr-l  + K 

- + 0(Ay*  ) 


2Ay* 


Equation  (21a) , 


o^1  - er  (er  + er,  - 29r) 

n n _ n-1  n+1  n 


At* 


Ay* 


iiS£.  mr_i  i + " 

St  dt*  1 n J S*r  dt* 


Vn_  dS*  (°n+I  ~ °n-l) 


2Ay* 


Equation  (21f), 

K-2  - ‘%r-l  + *>J> 


St 


2Ay* 


- a* 


s*m  - s*r 

At* 


0(  ) denotes  the  order  of  the  error  term  as  given  by  Ozisik,  (16). 
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